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ABSTRACT 



A general, modular, pulse propagation model for underwater acoustics that is based on 
linear systems theory for sound-speed profiles as a function of depth is presented. The 
development and computer implementation of the model, together with results from 
preliminary computer simulation studies involving the transmission of CW and LFM pulses 
from a planar array of complex weighted point sources is reported. The studies examined 
free-space propagation problems (i.e., no boundaries) in homogeneous and 
inhomogeneous media using a transfer function of the ocean medium based on the WKB 
approximation. The two main outputs from the model are the predicted complex acoustic 
field as a function of frequency and spatial location and the time-domain output electrical 
signal from each element in a receive planar array. 
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I. INTRODUCTION 



For small-amplitude acoustic signals it can be shown that the linear acoustic wave 
equation governs the propagation of sound in a fluid medium [Ref. l:pp. 23-31] [Ref. 
2:pp. 98-105]. Due to the linear nature of this equation, the ocean medium can be treated 
in general, as a linear, time-variant, space- variant, random filter or communication channel 
completely characterized by a transfer function [Ref. 3:pp. 1-6]. As a result, linear 
systems theory can be applied to the ocean medium propagation problem and through the 
use of coupling equations [Ref. 4] analytical expressions can be derived for the complex 
acoustic field and the output electrical signal at each element in a hydrophone array in terms 
of: 



• the frequency spectrum of the transmitted pulse 

• the far-field directivity function or beam pattern of the transmit array 

• the ocean medium transfer function 

• the far-field directivity function of the receive array. 

These expressions are developed for planar arrays, since they include linear arrays and 
single transducers as special cases, and are amenable to computer simulation studies. 

The application of Fourier analysis and linear systems concepts to problems in physics 
is not unique and has been applied extensively to the study of optics [Ref. 5]. What is 
unique is the application of these ideas to the physics of ocean acoustics and the treatment 
of the ocean medium, in general, as a time-variant, space-variant linear system [Ref. 3:pp. 
7-28] [Ref. 6]. 
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In this thesis, we will present a general, modular, pulse propagation model for 
underwater acoustics that is based on linear systems theory and the coupling equations for 
sound-speed profiles that are a function of depth. Since the coupling equations are the 
formal solution of the pulse propagation problem, and since the coupling equations depend 
on the transfer function of the ocean medium, the only thing that changes from problem to 
problem from an ocean acoustics point of view, is the ocean medium transfer function. 
Therefore, regardless of the problem under consideration, the coupling equations only need 
to be programmed once. Preliminary computer simulation studies can then carried out for 
the specific case of a free-space propagation problem using a transfer function based on the 
WKB approximation. 

The evaluation of the ocean medium transfer function is based on the full-wave solution 
technique, which involves the evaluation of wave propagation vector component integrals. 
This method is different from, but can be related to, ray acoustics and normal mode 
approaches to solving the propagation problem. As a consequence, our model of the pulse 
propagation problem is termed a full-wave solution. The computer implementation of a 
full-wave solution to the propagation problem is not new [Ref. 7] but an implementation of 
a model also based on linear systems concepts is. 

For the purposes of this thesis the transmit and receive arrays were assumed to be 
motionless which meant that we were able to treat the ocean medium as a linear, time- 
invariant, space-variant filter. Chapter 2 introduces and applies the linear systems 
approach to the time-invariant ocean medium problem and develops the coupling equations 
as far as possible without making any further simplifying assumptions about the ocean 
medium. As a consequence the results obtained can be applied to the propagation of an 
acoustic signal in any fluid medium. Chapter 3 then develops the coupling equations for 
the case of an ocean medium that is axially symmetric about the depth axis. This 
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development was based on the fact that axial symmetry was due to the speed of sound 
being a function of depth only, i.e., the case we wish to solve for. The results of this 
analysis can then be applied to any ocean medium description that exhibits axial symmetry. 
Chapter 4 then goes on to develop the ocean medium transfer function for the specific case 
of the ocean modelled with the WKB approximation. 

At this point, the development of the pulse propagation model is complete and is 
capable of producing two major outputs. The first output provides the magnitude and 
phase of the complex acoustic field incident upon the receive array as a function of 
frequency and spatial coordinates. This information is important, for example, for target 
localization using matched-field techniques. The second output is the time-domain output 
electrical signal at each element in the receive array. This information is important, for 
example, to illustrate pulse distortion due to the effects of dispersion in a waveguide, and 
as input to space-time signal processing algorithms. 

An integral pan of the model was the development, in Chapter 5, of a signal generator 
to simulate the transmitted electrical pulses that would be converted to acoustic energy for 
the propagation problem. It was assumed that the pulses were narrowband amplitude-and 
angle-modulated carriers. The signal generation scheme was based on a truncated Fourier 
series and complex envelope representation of narrowband signals. 

Chapters 6 and 7 contain computer simulation results for the specific case of the 
ocean medium that is characterized by the WKB approximation. These chapters deal with 
simulation results for an ocean medium characterized by a constant speed of sound and a 
sound-speed profile with a single constant gradient, respectively. 
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II. LINEAR SYSTEMS 



A. THE LINEAR SYSTEMS APPROACH 

Treating acoustic signals as small fluctuating disturbances in a fluid medium, it can be 
shown that the principles of fluid mechanics can be used to derive an equation governing 
the propagation of sound in a fluid medium. The linear acoustic wave equation, given by 

2 i 

V cp(t,r) - — cp(t,r) = x M (t,r), (2.1) 

c (r) 9t 

(where cp(t,r) is the velocity potential, in m2 / sec, at time t and spatial location r, x M (t,r) 
is the input acoustic signal to the medium, in 1 / sec, and c(r) is the speed of sound in the 
medium, in m / sec) can be derived by making a number of restrictive assumptions and 
applying them to the equations of state, continuity and momentum for a fluid. These 
restrictions allow us to derive a simple equation, of the form given by Eq. (2.1), which 
nonetheless is adequate to describe most commonly encountered acoustic phenomena. 
[Ref. l:pp. 23-31] [Ref. 2:pp. 98-105] [Ref. 3:pp. 1-3] 

The more commonly encountered acoustic quantities of acoustic pressure, p(t,r), in 
N / m2 and the acoustic fluid particle velocity, u(t,r), in m / sec can be obtained from the 
velocity potential as follows : 



p(t,r) = -p 0 ^cp(t,r) 



( 2 . 2 ) 



where p 0 is the constant equilibrium density of the fluid in kg / m3, and 
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u(t,r) = V(p(t,r). 



(2.3) 



It should be noted here that the concept of a scalar velocity potential arises from one of the 
restrictions used in the development of Eq. (2.1), namely that the particle motions 
associated with sound waves are irrotational. Therefore, the relationships given by Eqs. 
(2.2) and (2.3) are governed by this assumption. 

The linear acoustic wave equation given by Eq. (2.1) is so called because it is a linear 
partial differential equation relating the velocity potential and input signal. Since Eq. (2.1) 
can be thought of as a linear operation denoted by 



c (r) 3t 2 



(2.4) 



that relates the velocity potential and input signal, we can use a linear systems theory 
approach to describe the effect of the fluid medium on the propagation of an input signal 
through the medium. 

Using a linear systems approach, we describe the fluid medium as being a linear system 
that is characterized by an impulse response that is a function of both time and space. To 
then find the output (i.e., the acoustic velocity potential <p(t,r)) from such a system due to 
some input (i.e., the acoustic signal x M (t,r)), we simply convolve the input with the 
impulse response of the system. The basic geometry of this problem is illustrated in 
Fig. 2.1. 

This operation can also be done in the frequency ( in Hz ) and spatial frequency (in 
cycles / m) domains using the transfer function of the system. Since this is analogous to 
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( depth ) 

Fig. 2.1 Ocean Medium Geometry 

the concept of filtering, we say that our system is being modelled as a time-variant, space- 
variant, linear filter. [Ref. 3:pp. 3-6] 

Up until this point, we have been talking about any general fluid medium. We now 
wish to consider this problem with specific reference to the ocean medium. For the 
purposes of this thesis, we will model the ocean medium as a time-invariant, space-variant 
linear filter. Time invariance means that we will not be considering Doppler shifts between 
the input and output of the filter. In relation to the ocean medium we are modelling, this 
means that we are ignoring such things as motion between the transmitter and receiver, 
target motion and discrete point scatterers. 

Space variance means that there is a shift in spatial frequency between the transmitted 
and received signal. This means that sound speed is a function of position and there exists 
a non-uniform sound-speed profile. As a result, a transmit signal will scatter due to 
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refraction. If we think of a ray interpretation of acoustic propagation, the rays will bend 
rather than travel in straight lines as is the case for homogeneous media. 

B. THE LINEAR, TIME-INVARIANT, SPACE- VARIANT FILTER 

Consider a system modelled by a linear, time-invariant, space-variant filter. The time- 
invariant, space-variant impulse response of the filter, h(x,r 0 ;r), will then completely 
characterize the system being modeled. The relationship between the input x(t,r) and the 
output y(t,r) of such a filter is given by 



y(t.r) = 




x(t-x,r-r 0 )h(x,r 0 ;r)dx dr 0 



(2.5) 



where the function h(x,rQ;r) is the response of the filter positioned at spatial location 
r = (x,y,z), due to an impulse applied at position r 0 = (x 0 ,y 0 ,z 0 ), x seconds ago. For 
the specific case of modeling the ocean medium, the terms y(t,r) and x(t,r) in Eq. (2.5) 
would equate to <p(t,r) and x M (t,r) in Fig. 2.1, respectively. 

As an input to our filter, consider a signal x(t,r) of the following form: 



x(t,r) = e- 



j27ift -j2;rvr 



( 2 . 6 ) 



where f is the input frequency in Hz and v = (f x ,f y ,f z ) ^ i n P ut spatial frequency 
components in cycles / meter. The spatial frequencies are a function of direction and 
wavelength given by 



f x =uM , 



(2.7a) 
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fy = V/X 



(2.7b) 



and 



f z =w/X (2.7c) 

where u, v and w are the direction cosines with respect to the positive X, Y and Z axes. 
The surfaces of constant phase, commonly referred to as wavefronts, for the signal 
described by Eq. (2.6) are given by 27tv*r = constant, i.e. a plane., Consequently, Eq. 
(2.6) describes a time-harmonic plane wave traveling in the direction of the propagation 
vector k = 27tv, normal to the wavefront. [Ref. 2: pp. 107-108] 

To find the filter output due to a time-harmonic plane wave input, we combine 
Eqs. (2.5) and (2.6) which yields 

y(t,r) = e j2nfl e^ 2rcv ’ r H(f,v;r), (2.8) 



where 



H(f,v;r) 



■a: 



h(T,r 0 ;r)e- j2,rfx e j27CV * r odT 



dr n 



(2.9) 



is the transfer function of the linear, time-invariant, space-variant system (or filter). It is 
completely analogous to the transfer function H(f) of linear time-invariant (LTI) systems 
often encountered in electrical circuit theory and signal analysis. 
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Equation (2.9) is in the form of a multidimensional Fourier transform of the impulse 
response and can be written as 



H(f,v;r) = F t F r { h(x,r 0 ;r) } 



( 2 . 10 ) 



where F x represents the forward time-domain Fourier transform given by 



H(f,r 0 ;r) = F T { h(x,r 0 



-/■ 



h(x,r 0 ;r)e' j2,rft dx , (2.11) 



and F rQ represents the forward spatial Fourier transform given by 



■/: 



H(x,v;r) = F r ( h(T.r 0 ;r) ) = | h(T,r 0 ;r) e j2,tv-r ° dr 0 . (2.12) 



Comparing Eqs. (2.8) and (2.9), we see that we have two methods of obtaining the 
transfer function of our filter. Since the transfer function completely characterizes the 
behavior of the filter and, hence, the ocean medium that the filter is modeling, we need 
some method of determining this quantity. Equation (2.9) is not much use for our 
problem since it assumes that we already know the impulse response of the filter which we 
do not. However, Eq. (2.8) is much more useful. It says that if we apply a plane wave 
input to the system, then there is a simple relationship between the input and output that 
determines the value of the transfer function. However, this relationship is more versatile 
than may first appear. It tells us that if we can express an output due to any type of input 
as the product of the time-harmonic plane wave given by Eq. (2.6) and some other factor, 
then this other factor is the transfer function. 
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Expressing the impulse response term in Eq. (2.5) in terms of the two dimensional 
inverse Fourier transform of its transfer function and combining with the input term, it can 
be shown that the two dimensional Fourier transform of the output of a linear time- 
invariant, space-variant filter, given by Eq. (2.5), has the form 



which is the frequency and angular spectrum of the output. Note that the output spatial 
frequency term is different, in general, from that at the input, and that Y * XH . The 
exception to this rule is when the ocean medium is homogeneous. In this case the filter 
model is space-invariant and Eq. (2.13) collapses to the form Y = XH , which is the same 
as the familiar result encountered for one-dimensional LTI systems. 

We can determine the impulse response once we know the transfer fountain by using 
the inverse transform given by 




(2.13) 



h(x,r 0 ;r) = 




(2.14) 



which can be expressed in transform notation as 



h(t,r 0 ;r) = F f 1 Fv { H(f,v;r) } 



(2.15) 



where Ff 1 represents the inverse time-domain Fourier transform and Fy represents the 
inverse spatial Fourier transform. The form of Ff 1 and Fy is the same as for Eqs. (2. 1 1) 



and (2.12) with the sign of the exponent changed and the variable of integration changed to 



f and v, respectively. This is the same as for the familiar one-dimensioned time-frequency 
Fourier transform. 

Note that the development in this section has so far been completely general and can be 
applied to any linear, time-invariant, space-variant system. We have not made any 
assumptions regarding the space variance of the system, i.e., we have not specified a 
particular form for the sound speed-profile as a function of three-dimensional space. 

C. THE PROBLEM STATEMENT 

Let us now consider the underwater acoustic problem from the viewpoint of linear 
systems theory. Our aim is to be able to compute the output electrical signal from a receive 
transducer due to an applied electrical signal at a transmit transducer. Simply stated, this 
process consists of three steps : 

• conversion of an electrical signal into a transmitted acoustic signal, 

• propagation of an acoustic signal through the ocean medium, and 

• conversion of a received acoustic signal into an electrical signal 

Our intention is to model each of these steps as a time-invariant, space-variant / invariant 
filter as appropriate. The system model conforming to this plan is shown in Fig. 2.2 in 
block diagram form. 

The remainder of this chapter will deal with the development of expressions describing 
the behavior of the transmit and receive arrays, and the relationship or coupling equation 
that relates the output, y(t,r), to the input, x(t,r). The development of the transfer function 
that characterizes the ocean medium will be discussed in the next chapter. 
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Fig. 2.2 System Model 
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D. COUPLING THE TRANSMITTED SIGNAL TO THE MEDIUM 

Consider now the specific problem of transmitting an electrical signal through the ocean 
medium. The electrical signal will be used to drive some transducer (or array of 
transducers) which will convert the electrical energy into acoustic energy. We therefore 
need to find a relationship that describes this transformation of an electrical signal x(t,r) 
into an acoustic signal x M (t,r) . The acoustic signal represents the rate at which fluid 
volume is added at time t and position r per unit volume of fluid. 

Suppose we model an infinitesimal region of some arbitrarily shaped transducer as a 
linear filter with a corresponding impulse response. The output from such a transducer at 
a given position r would be given by the time-domain convolution of the corresponding 
impulse response with the input signal. Translated into the frequency domain via a Fourier 
transform with respect to time, this relationship can be expressed as 

X M (f,r ) = X(f,r ) A T (f,r) , (2.16) 

where A T (f,r) is the complex frequency response (or aperture function) of the transducer 
at spatial location r. Taking the spatial Fourier transform of both sides of Eq. (2.16), it 
can be shown that : 
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X(f,a)D T (f,v-a) da 



(2.17) 




where D t (f,a) is the spatial Fourier transform of the complex frequency response of the 
transducer and is known as the far-field directivity function or beam pattern. Note that the 
result given by Eq. (2.17) is valid for a completely arbitrary volume aperture. 

If we make the simplifying assumption that an identical electrical signal is applied at all 
spatial locations, r, of the transducer then 



Taking both the time and spatial domain Fourier transforms of a signal of the form of 
Eq.(2.18) and substituting into Eq. (2.17) yields 



If the transmit aperture is a planar array, then for most practical situations we can consider 
that an identical input electrical signal is applied to all the elements of the array before any 
electronics used for beamsteering. As a result, the restriction imposed by Eq. (2.18) does 
not impose any limitations. What we gain, however, when comparing Eqs. (2.17) and 
(2.19), is a significant simplification in our analysis. [Ref. 4:p. 4] 

E. THE TRANSMIT ARRAY 

Assume that the transmit aperture is a planar array of identical transducer elements. 
The product theorem for planar arrays, given by 



x(t,r) = x(t). 



(2.18) 



X M (f,v) = X(f)D-j<f,v). 



(2.19) 
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D(f,v) = S(v)E(f,v), 



( 2 . 20 ) 



states that the far-field directivity function of such an array is equal to the far-field 
directivity function of one of the transducer elements, E(f,v) , multiplied by the far-field 
directivity function of an equivalent array of complex weighted point sources, S(v) . We 
will assume for the purposes of this thesis that our transmit array is located in the XY plane 
and consists of MT x NT (odd) complex-weighted point sources. We will further assume 
that the array elements are equally spaced in the X and Y directions; however, it should be 
noted that this is not a constraint imposed by the use of the product theorem. 

Since the spatial Fourier transforms deal with position relative to some source, it 
makes sense to locate the transmit array in a convenient position to simplify the analysis. 
Accordingly, we position the array such that its center forms the origin of the coordinate 
frame of reference in the X and Z directions while the ocean surface forms the center in the 
Y direction, i.e., the array is centered at (x 0 = 0,y o = yx,z 0 = 0) . The geometrical 
interpretation of this configuration is shown in Fig. 2.3, which is also the physical 
representation of the block diagram of Fig. 2.2. 

Given these conditions, Ziomek [Ref. 3:p. 126] shows that the far-field directivity 
function of the transmit array is given by 

D T (f,V)= £ X (22I) 

m = -MT n = -NT 



where c mn (f) is the frequency dependent complex weight at the transmit element (m,n), 
and d X T and d yr are the interelement spacings in the X and Y directions, respectively. 
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( depth ) 



Fig. 2.3 Transmit and Receive Array Geometry 



The last complex exponential is a phase factor that accounts for the array not being centered 
at the origin. Also, 



MT = ( MT - 1 ) / 2 



( 2.22a) 



and 



NT = ( NT - 1 ) / 2. (2.22b) 

It should be noted from Eq. (2.21) that a single omnidirectional point source and a linear 
transmit array are special cases of the planar transmit array. This is also what we would 
expect from examining this problem from a purely physical point of view. If we assume 
that the transmitted electrical signal is applied to the array elements before the complex 
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weighting, then we can say that an identical electrical signal is applied to all the elements in 
the array. Accordingly, Eq. (2.21) can be combined with Eq. (2.19) to describe the 
acoustic input to the ocean medium. 

Finally, what is meant by the far-field directivity function? It can be shown that a field 
point at a range r is considered to be in the far field if 




(2.23) 



where R is the maximum radial extent of the aperture. For our planar transmit array, this 
distance is given by 

R = [(Mrd XT ) 2 +(Nrd YT ) 2 l . (2.24) 

The physical significance of being in the far field is that the range is large compared to a 
wavelength and the separation between the transducer elements. Consequently, any 
second or higher order terms in the expression for the directivity function can be ignored 
and the mathematics of the problem is simplified. [Ref. 3:p. 38] [Ref. 2:p. 170] 

F. THE PHASED ARRAY 

Now let us take a closer look at the frequency-dependent complex weighting function 
Cmn (0 of the far-field directivity function given by Eq. (2.21). Since this weighting 
function is complex, it can be expressed as 

c mn (0 = a mn (0 mn (2.25) 
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which consists of an amplitude weighting component a mn (f) and a phase weighting term 
®mn (0 . Now consider that this phase weighting is a linear phase variation of the form 

®mn(f) = ( fx^XT + fy^YT) (2.26) 

where m, n, d xx , and dyr are defined as for Eq. (2.21). If we define the spatial 
frequency terms in Eq. (2.26) as 



fx=u'A (2.27a) 

and 

f Y =v'M, (2.27b) 

then, substitute Eqs. (2.25) and (2.26) into Eq. (2.21), it can be shown that the resulting 
far-field directivity function has been steered in the direction u = u’ and v = v'.[Ref. 3:pp. 
132-134]. 

The term phased array means nothing more than an array which employs beamsteering 
using phase weighting. Inspection of our assumed form of the phase weighting term 
shows that it is both a separable function and an odd function of the indices m and n. 
These are necessary conditions in order to carry out the beamsteering. It should be noted 
from the preceding development that no assumption was made as to whether the amplitude 
term was either separable or an odd function of the indices m and n. For the purposes of 
this thesis, we will only be considering rectangular amplitude weighting, which is a 
separable function and leads to a simplification in the analysis. 
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Since the phase weighting is already a separable function, the assumption that the 
amplitude weighting is also separable results in the entire complex weighting term given by 
Eq. (2.25) being separable. Substitution of this form of weighting into Eq. (2.21) results 
in the expression for the far-field directivity function being separable which gives 



Dj(f,v) = Dy(f,f x ) Dy^fy) 



(2.28) 



where 



MT’ 

D T (f,fx)= I Ae iM x mdxT e- j 2 ,rf x md XT i 

m=-MT' 



(2.29) 



A is the magnitude of the amplitude weighting, and a similar expression exists for 
D T (f,f Y ). To give a more intuitive feel for what this far-field directivity function looks 
like, we can express Eq. (2.29) in the alternative form 



sin[7t(fx - fY)MTdxT] 
d t <f,f x ) = a x x ; — — 

sin[7t(f x -f x ) d XT] 



(2.30) 



which looks similar to a sine function when viewed in direction cosine space. The 
expression for D T (f,f Y ) is of the same form as Eq. (2.30) except that it is multiplied by the 
last complex exponential phase factor in Eq. (2.21) that accounts for the array not being 
centered at the origin. 
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G. RECEIVED ACOUSTIC FIELD 



We next need to calculate the acoustic signal at some spatial location r = (x,y,z) after 
transmission through the ocean medium. Applying the space-variant, time-invariant 
impulse response of the ocean medium to Eq. (2.5) yields 



If we express the input signal in the previous expression in terms of its frequency and 
angular spectrum and collect terms which combine with the impulse response to produce 
the transfer function of the ocean medium, we obtain 



which is the frequency spectrum of the output acoustic field and is one of our desired 
results. 

Next we need to consider some receive aperture to convert the acoustic energy into 
electrical energy. For the purposes of this thesis, we will assume that the receive aperture, 
like the transmitter, is a planar array that is characterized by a complex aperture function 
A R (f,r). Analogous to the development of Eq.(2.16), it can be shown that the 
corresponding output electrical signal is given by 



y M (t,r)= x M (t-T,r-r 0 )h M (x,r 0 ;r)dTdr 0 . (2.31) 





(2.32) 




(2.33) 
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We will assume that our receive array is located in a plane parallel to the XY plane and 
consists of MR x NR ( odd ) point receivers centered at (x R ,y R ,y R ). We will further 
assume that the array elements are equally spaced in the X and Y directions. Due to the 
assumption of point receivers, the aperture function can be expressed as 

MR' NR' 

A R (f,r)= X X w mn (f)8[x-(x R +md XR )]8[y-(y R +ndYR)]5[z-z R ] (2.34) 
m=-MR' n=-NR' 

where w mn (f) is the frequency dependent complex weight at element (m,n), dxR and dyR 
are the interelement spacings in the X and Y directions, respectively, 

MR’ = ( MR - 1 ) / 2, (2.35a) 



and 



NR' = ( NR - 1 ) / 2. (2.35b) 

For this thesis, we are not concerned with processing the received electrical signals. 
What we want to accomplish is the ability to compute the complex acoustic field at the 
receive array as a function of frequency and spatial coordinates, Eq. (2.32), and the time- 
domain electrical signal at each element in the receive array. The signal processing that 
will be done on the computer simulated received electrical signals will be via auxiliary 
programs, e.g., FFT and adaptive beamformers. Accordingly, we are primarily concerned 
with formatting our outputs in a form acceptable to those programs. 
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As a consequence of these requirements, and in an attempt to minimize computation 
time, we will assume the simplest form of the aperture function where the complex 
weighting function w mn (f) is unity. Applying this simplifying assumption to Eq. (2.34) 
and substituting the resulting form of the aperture function into Eq. (2.33) yields 



y(t,m,n,z R ) = 




Y M (f,m,n,z R )e i2nfl df 



(2.36) 



which is the time domain electrical signal at element (m,n) and is our other desired result 

H. OVERALL SYSTEM COMPLEX FREQUENCY RESPONSE 
If we combine Eqs. (2.19) and (2.32) we obtain 

Y m (f,r) = X(f)H(f,r) (2.37) 



where 



H(f,r) = 




D T (f,v)H M (f,v;r)e' j27tv ’ r dv 



(2.38) 



is known as the overall system complex frequency response [Ref. 4:p. 5] [Ref. 6:p. 1672]. 
Note that in the development carried out in this chapter we have not made any assumptions 
about the ocean medium transfer function. The only assumption used in the derivation of 
Eq. (2.35) was that an identical electrical signal was applied at all spatial locations of the 
transmit array before the application of the complex weights. 
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Taking the inverse time-domain Fourier transform of Eq. (2.37), and using the 
expansion of Eq. (2.38) it can be shown that 

y M (t.r)=J J X(f)D T (f,v)H M (f,v;r)e‘ j2,K ' r e i2,tft dvdf (2.39) 



which is the coupling equation that relates the output acoustic signal from the medium, 
Ym ( t ’ r ) , to the frequency spectrum of the applied elctrical signal X(f). 
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III. THE AXISYMMETRIC CASE 



A. OVERALL SYSTEM COMPLEX FREQUENCY RESPONSE REVISITED 
Let us examine the expression for the overall system complex frequency response, 
given by Eq. (2.38), that was developed in the previous chapter. It would appear from an 
initial inspection that Eq. (2.38) is a three-dimensional integration with respect to 
dv = df x df y dfz . However, it should be noted that since the integrand is a function of 
the frequency, f, as well as the spatial frequencies, v = ( fx>fY>fz)> the integral can be 
reduced from a three-dimensional to a two-dimensional integral. This is due to the manner 
in which the spatial frequencies given by Eq. (2.7) are related to the frequency f via 
direction cosines. The property of direction cosines that the sum of their squares always 
equals unity means that, since we know f, we can express one of the direction cosines in 
terms of the other two. Thus, we do not need to integrate with respect to all three spatial 
frequencies in order to span direction cosine space. Using the spatial frequency definitions 
of Eq. (2.7), it can be shown that 

xe -j2,( fx x + f Y y + f z z )dfxdfz (3i) 



H(f,x,y,z) 



=/:/: 



where 



f Y =±( (f/c 0 ) 2 -(fx + fl) ] , (fx +4) fi(f/c„) 2 , (3.2a) 



and 



23 



(3.2b) 



f Y = +il(fx+fz)-(f/Co) 2 l * < (fx +fl)> (f/c 0 ) 2 , 

and c 0 is the speed of sound at a particular transmit array point source element. The plus 
(minus) sign in Eq. (3.2a) is chosen whenever y > yo (y < yo) and corresponds to the 
field point being deeper (shallower) than the source point, i.e., downward (upward) 
traveling waves. The minus (plus) sign in Eq. (3.2b) corresponds to the plus (minus) 
sign of Eq. (3.2a) in order to generate evanescent waves. 

Note that we have not made any assumptions about the ocean medium transfer function 
in our development up until this point. However, for the purposes of this thesis, we are 
going to assume that the acoustic energy of the transmitted signal is contained within a 
sound channel, e.g., the SOFAR or deep sound channel. In this situation the sound 
propagates between upper and lower turning points with no surface and bottom boundary 
interaction. We therefore expect the acoustic field trapped within this channel to exhibit 
some form of cylindrical spreading. If this assumption is correct, then cylindrical 
coordinates should provide an easy reference frame against which to describe this behavior. 

B. THE CYLINDRICAL COORDINATE SYSTEM 

The cylindrical coordinate system used in this thesis is shown in Fig. 3.1, where the 
transformation equations that relate the rectangular and cylindrical coordinates are given by 

x = rcos<]), y = y, and z = r sin<]) (3.3) 

2 2 ~ 

x + z is the polar radial distance. 

In addition, the spatial frequencies must be transformed using the following set of 
equations ( see Fig. 3.2 ): 
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z 




Y 

( depth ) 



Fig. 3.1 The Cylindrical Coordinate System. 



2nf z 
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f x = f r cos£, fy =fy, andfz = f r sin£ 



(3.4) 



V 2 2 

f x + f z is the spatial frequency in the radial direction. As a result. 



df x df z ->f r df r dC. 



(3.5) 



The physical significance of the two different angles in Eqs. (3.3) and (3.5) is that <j> relates 
to the position of some field point whereas £ relates to the direction of propagation, and 
they are not necessarily the same. 

Applying these definitions to Eq. (3.1), it can be shown that the overall system 
complex frequency response can be expressed in cylindrical coordinates as 




x e "j27t(f r r cos£ sin<> + f Y y + f r r sin£ sin<(>) ^ ^ 



where 




(3.7a) 



and 




(3.7b) 



and Co is the speed of sound as previously defined. 
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Now consider the specific transmit array geometry discussed in Chapter 2. If we 
substitute the far-field directivity function for this array, as given by Eq. (2.21), into the 
previous equation, and perform the necessary transformation from rectangular to 
cylindrical coordinates, we find that 



x e -j2*f r (r cos; sin0 + r sin; sin<j> + md xx cos;) e -j 27 tf Y (y - y T - ndyiOf df d £ (3 ^ 

C. THE AXISYMMETRIC OCEAN MEDIUM TRANSFER FUNCTION 

In this thesis, we shall assume that the speed of sound in the ocean medium is a 
function of depth only. It will be shown in the following chapter that under this condition 
the ocean medium transfer function is axi symmetric, i.e., independent of the azimuthal 
angle ; and, as a result, it can be expressed in the form 



where y 0 =yj - ndyj is the depth of a point source in the nth horizontal row of the 
transmit array. Note the presence of the terms y and y Q in Eq. (3.9). These terms 
indicate that the ocean medium transfer function is dependent on the speed of sound at the 
position of the transmitter as well as the receiver. Since the assumed form of our 
transmitter is a planar array, we need to consider the depth and, hence, speed of sound at 
each element in the array when calculating the acoustic field at some spatial location r. This 
makes physical sense when considering Eq. (3.2), which shows that the speed of sound at 
the transmit element affects the direction of propagation and the threshold between 




H M (f,f r ,C,f Y ;r,<)),y) = H M (f,f r ,y 0 ;y) 



(3.9) 
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propagating and evanescent waves. If the form of the ocean medium transfer function 
given by Eq. (3.9) is substituted into Eq. (3.8), then the integrals can be regrouped to give 



oo MT' NT 



-j27tf Y (y - y 0 ) f 

' M 



HOD- I 1 I I 



c mn( f )H M (f,f r ,y 0 ;y)e' 



r 



0 m = -MT n = -NT 




-j27tf r (r cos£ sin<j) + r sinC sin<J> + mdxTCOsQ 



d£df r . (3.10) 



Next, consider the integration with respect to C in the preceding equation. Due to 
our placement of the transmit array in the coordinate system, the angle <j> is the horizontal 
polar angle (relative to the positive X axis) between the center of the transmit array and the 
field point of interest. Now define a new angle <|> m that is the horizontal polar angle 
between a specific array element (m,n) and the field point, as shown in Fig. 3.3. Note that 
the field point in Fig. 3.3 is in three-dimensional space, but that it is only the (x,z) 
component of that position that determines <j). Also note that since the angle is in the XZ 
plane, it is independent of the vertical position y and, hence, is independent of the array 
index n. Expressing the horizontal polar angle <}) in this integral in terms of the new angle 
<t> m and collecting terms, results in 

f 2n e -J 2rrf r( r cosC sin<j> + r sin£ sin<|> + md XT cosQ d £ f 2 * e -j 2jrf r r m c0S (C - <t> m ) d £ (3 j ^ 

*'0 •'0 



where 



r m =[(x-mdxT) 2 +z 2 ] 



1/2 



(3.12) 
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(m,n) 




Z 



(x,y,z) 



X 



Fig 3.3 Graphical Interpretation of <> m ( Plan View ). 



is the polar radial distance between element (m,n) in the transmit array and the field point, 
and mdxr is the horizontal distance between element (m,n) and the center of the array. 

Brekhovskikh [Ref. 8:p. 76] and Officer [Ref. 9:p. 126] show that the right-hand side 
of Eq. (3.1 1) is in the form of the integral representation of the zero-order Bessel function 
of the first kind where 




-j27tf r r m cos(C - <M 



'dC = 2jcJ 0 (2jtf r r m ). 



(3.13) 



Substituting Eqs. (3.1 1) and (3.13) into Eq. (3.10) yields 



00 MT' NT 



H(f,r) = 27t j X S c mn (0H M (f,f r ,y o ;y) J 0 (27tf r r m ) 



0 m = -MT n = -NT 



m = -MT n = -NT 




(3.14) 
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which is the fundamental result of this chapter. Comparing the previous equation with Eq. 
(3. 1), it can be seen that our expression for the overall system complex frequency response 
has been reduced from a double to a single integral. This was based on the assumption 
that the speed of sound was a function of depth only and, as a result, the corresponding 
ocean medium transfer function was axisymmetric. Since numerical integration is a time 
consuming process, it is surmised that the reduction in the number of integrations required 
in the computation of our computer-generated acoustic output should result in the saving of 
a significant amount of computer time. 

It should be emphasized that Eq. (3.14) is valid for any axisymmetric ocean medium 
transfer function. Since this condition is satisfied for the case when the sound speed is a 
function of depth only, Eq. (3.14) represents the general form of the overall system 
complex frequency response that we will use in the development of the next chapter. 
Remember that we are trying to develop a general, modular, pulse-propagation model 
based on sound-speed profiles that are a function of depth. Thus Eq. (3.14), when 
combined with the coupling equations of Chapter 2, is the general solution for this 
assumed type of sound-speed profile. In the next chapter we will develop a specific 
solution based on the WKB approximation. 

D. AXIAL SYMMETRY OF THE TRANSMIT ARRAY 

The only assumption made regarding axial symmetry in this development was that the 
ocean medium was axisymmetric. This makes physical sense when considering the ocean 
medium itself, in the absence of any source of acoustic signal. Since the speed of sound 
was assumed to be a function of depth only, there is nothing in the medium that should 
cause any azimuthal <)> dependence on the propagation of an acoustic signal. Accordingly, 
the ocean medium transfer function should be independent of the azimuthal angle <|). 
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The assumption of axial symmetry may at first glance appear to be invalid when we 
include the transmit array in the problem, as the beam pattern or far-field directivity 
function of the planar array is not generally symmetrical about the Y axis. However, 
consider the problem from the following perspective. From linear systems theory, the far- 
field directivity function of the array is the linear superposition of the far-field directivity 
functions of the individual array elements. The far-field directivity function of a point 
source is unity, i.e., it is omnidirectional. The far-field directivity function of a complex 
weighted point source is also omnidirectional, since the complex weight factor is merely a 
scaling term. Equations (2.20), (2.21) and (2.29) demonstrate that the far-field directivity 
function of the transmit array is a result of the way these omnidirectional beam patterns are 
combined rather than by altering the beam pattern of the individual point sources. The 
effect of the complex weighting is to control the way these individual far-field directivity 
functions are combined. 

Examination of Eq. (2.21) also shows that there is an additional phase term that is 
associated with each array element that accounts for that element's displacement from the 
origin of the coordinate axis system. Thus, we can assume that the far-field directivity 
function of the transmit array is the linear combination of complex weighted axisymmetric 
point sources. Accordingly, the overall system complex frequency response solution 
given by Eq. (3.14) can be thought of as applying to each individual array element, with 
the solution for the array being a linear combination of these results. This is nothing more 
than rewriting Eq. (3.14) with the summations outside of the integral. Since the 
integration and summations are with respect to different quantities, regardless of whether 
the integration is inside the summations or outside, the value of the overall system complex 
frequency response will be the same. 
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IV. THE OCEAN MEDIUM TRANSFER FUNCTION 



A. THE HELMHOLTZ WAVE EQUATION 

Examining the various coupling equations developed in Chapter 2, we see that they 
already incorporate an input acoustic signal x M (t,r), with corresponding frequency and 
angular spectrum X M (f,v). This means that the transfer function that characterizes the 
ocean medium, H M (f,v;r), is independent of the input as we would expect from linear 
systems theory. Accordingly, when trying to find a transfer function to describe the 
ocean medium we only need to solve the wave equation given by 



instead of the linear acoustic wave equation given by Eq. (2.1). 

The first assumption that we will make in finding a solution to Eq. (4.1) is that the 
velocity potential, <p(t,r), has a time-harmonic dependence, that is. 



If this form of the velocity potential is then substituted into Eq. (4.1) and the time-harmonic 
dependence is factored out, the result is the time-independent Helmholtz wave equation 
given by 




(4.1) 



tp(t,r) = tp(r)e j27lft 



(4.2) 



V 2 <p(r)+k 2 (r)<p(r) = 0 



(43) 
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where k(r) = 27tf /c(r) is the wave number and c(r) is the speed of sound in the fluid 
medium at spatial location r. 

As discussed in the previous chapter, we will be assuming for the purposes of this 
thesis that the speed of sound is a function of only one spatial variable, the depth y, i.e., 
c(r) = c(y). While this assumption is not in general true, it is nevertheless adequate 
enough to describe the ocean medium for most circumstances where we are not considering 
propagation ranges longer than say, 20 km, or propagation across significant 
oceanographic perturbations such as fronts or eddies. We can also consider that the 
sound-speed profile in the region of the deep sound channel is less variable than for 
shallow sound channels or surface ducts. Thus, the assumption that the sound-speed 
profile is not dependent on geographical position is adequate for our model, within the 
limitations discussed. 

This assumption allows us to make the substitution k(r) = k(y) in Eq. (4.3) which 
greatly simplifies the analysis. Since k(r) is now a function of only one variable, we can 
use the method of separation of variables to solve for the acoustic velocity potential. 
Applying this solution technique in cylindrical coordinates, we assume a solution for cp(r) 
of the form 



<p(r) = R(r)<D«}>)Y(y) (4.4) 

where r = (r,<J>,y) is the spatial location in cylindrical coordinates (see Fig. 3.1) and the 
Laplacian appearing in Eqs. (4.1) and (4.3) is defined as being 

,4.5, 

dr 2 r dr r 0<j> dy 
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B. AZIMUTHAL ANGLE AND RANGE TERM SOLUTION 



Since the speed of sound is a function of depth only, there is nothing in the medium 
that should cause any azimuthal, i.e., <J>, dependence on the propagation of an acoustic 
signal. Accordingly, the solution for the acoustic velocity potential should be independent 
of the azimuthal angle <j>. This is referred to as the ocean medium having axial symmetry. 
Applying this assumption, Eq. (4.4) reduces to 

cp(r) = R(r)Y(y) (4.6) 



and the partial differential with respect to azimuthal angle in the Laplacian of Eq. (4.5) 
drops out. 

If we substitute Eq. (4.4) into Eq. (4.3) and assume that k(r) = k(y), the method of 
separation of variables yields 



^M + iM ttrR(r)=0 

dr r dr 



(4.7) 



where the wave number in the radial direction, k r , has been used as a separation constant. 
Making the substitution R(r) = g( k r r ) into the previous equation, it can be shown that Eq. 
(4.7) reduces to Bessel's differential equation of order zero. The solution to this equation 
can be expressed in either of the forms 

R(r) = A J 0 (k r r) + B Y 0 (k r r) (4.8a) 



or 
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(4.8b) 



R(r) = AHf, 1) (k I r) + BHS, 2) (k r r) 



where Jq and Yq are the zero order Bessel functions of the first and second kind, 
respectively, with Yq also being known as the zero order Neumann function. The 



respectively, or Bessel functions of the third kind. Since Eqs. (4.8a) and (4.8b) are 
alternative expressions describing the same quantity, the Hankel functions should be able to 
be expressed in terms of the Bessel functions of the first and second kind and visa-versa. 
This is, in fact, the case and we will make use of one of these relationships at a later stage 
of the analysis. 

Since there are two alternative representations of the range term R(r), we would expect 
that one representation is more useful than the other for the situation that we wish to model. 
If we consider the approximations for large argument for the two Hankel functions 



then because of our assumed form of the time-harmonic dependence of the velocity 
potential given by Eq. (4.2), Eq. (4.9a) represents an incoming wave while Eq. (4.9b) 
represents an outgoing wave. Since we are assuming that our transmitted acoustic signal 
is traveling in a channel with no boundary interaction, we will have only outgoing waves 



functions Hq" and Hq^ are zero order Hankel functions of the first and second kind. 




(4.9a) 



and 




(4.9b) 
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from the source and no incoming waves due to reflection. It then appears that Eq. (4.8b) 
with the constant A set to zero, is an appropriate solution for our range term. We will 
therefore assume that our range term solution is of the form 

R(r) =B r H^ ) (k r r) (4.10) 

where B r is a constant. By comparison, if we look at the approximations for large 
argument of the Bessel functions of the first and second kind, we find that they are 
comprised of both incoming and outgoing wave components and so are not well suited to 
our problem. 

Upon substitution of Eq. (4.10) into Eq. (4.7) however, we find that our assumed 

form of the solution for the range term is a solution everywhere except at the origin, i.e., 
( 2 ) 

r = 0, since Hq "blows up" at this point and equals -j ~ . Using a volume integral 
approach, Gauss's theorem, and the asymptotic expression for the Hankel function given 
by 



lim Ho 2) (k r r) =— In (k r r), (4.11) 

r — » 0 7lr 

it can be shown that Eq. (4.10) is a solution to 

[V 2 +k r 2 ]R(r) = 2^£.8<r) (4.12) 

for all r, where Eq. (4.12) is the same as Eq. (4.7) with the addition of a forcing function 
or source term. This makes physical sense when we remember that Eq. (4.10) describes 
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outgoing waves in the radial direction. Since these outgoing waves have to come from 
somewhere and there are no incoming waves, there must be a source of some description, 
that generates the waves. Thus Eq. (4.10) is a solution to Eq. (4.12) that satisfies the 
boundary condition of a point source at the origin. The choice of a point source as the 
forcing function is appropriate since we have assumed that our transmit aperture is a planar 
array of point sources. 

C. THE WKB APPROXIMATION 

For the depth dependent term Y(y) of our solution for the velocity potential, we will 
obtain an approximate solution based on the geometrical optics approximation, otherwise 
known as the WKB (Wentzel, Kramers, and Brillouin) approximation. If we substitute 
Eq. (4.4) into Eq. (4.3) and assume that k(r) = k(y), the method of separation of variables 
yields 



where the wave number in the vertical (or depth) direction, ky (y), is calculated using the 
relationship 



If the vertical wave number was a constant, i.e., ky(y) = ky, meaning that it was no 
longer a function of depth, then Eq. (4.13) would have the exact solution 




(4.13) 



ky(y) = k 2 (y)-k r 2 . 



(4.14) 



Y(y) = Ae- jkYy +Be jkYy 



(4.15) 
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If we now assume that the local medium is homogeneous, that is the variations of the 
properties of the medium per wavelength are small, then a solution similar in form to Eq. 
(4.15) would be a good approximation. The WKB method assumes that the solution to 
Eq. (4.13) is of the form 



Y(y) = A Y (y)e jeY(y) (4.16) 

where A Y (y) and 0 Y (y) are real amplitude and phase functions respectively. Substituting 
this assumed form of our solution into Eq. (4.13), and making the simplifying assumption 
that 



I A Y "(y) / A Y (y) | « [© Y '(y)] 2 



(4.17) 



yields 



Y(y)~ 



V|k Y (y)| 



r -j l 


f y k Y (0dC jf 


"* k Y (0<X 




[ A Y e J 


y D + B Y e J ■ 


J 


, (4.18) 



which is the approximate WKB solution [Ref.3:pp. 208-220] [Ref. 8:pp. 129-134]. It 
should be noted that the only assumption used that leads to this not being an exact solution 
is that given by Eq. (4.17). This method also assumes that there is no reflection in a 
horizontally stratified medium. 

The LHS of Eq. (4.17) is the relative differential of the time rate of change of 
amplitude and is a measure of the rate of change in the properties of the medium. If the 
medium was homogeneous, then this term would be zero since the characteristics of the 
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medium would not be varying. Therefore this term is a measure of the inhomogeneity of 
the medium. The RHS of Eq. (4.17) is inversely proportional to wavelength since for a 
given frequency the rate of change of phase increases with decreasing wavelength. 
Equation (4.17) is therefore a reasonable simplifying assumption since it is nothing more 
than a mathematical formulation of the assumption made earlier that the variations of the 
properties of the medium are small compared to a wavelength. 

The integrals in Eq. (4.18) give the phase change as the waves propagate between 
depths y and y<> while the factor in front is to preserve the conservation of energy. Due to 
the assumed form of the time-harmonic dependence of the velocity potential given in Eq. 
(4.2), the first and second terms in Eq. (4.18) correspond to waves traveling in the positive 
(down) and negative (up) Y directions, respectively. It should be noted that the WKB 
solution is then the superposition of two waves traveling without interaction in opposite 
directions. Also note that the approximate solution is not valid as k y (y) approaches zero 
since the solution approaches infinity. The point at which this occurs is referred to as a 
turning point. 

D. THE VELOCITY POTENTIAL SOLUTION 

Combining Eqs. (4.6), (4.10) and (4.18), we obtain a solution for the velocity potential 
of the form 




(4.19) 



where 




(4.20a) 
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and 




which is our separable function solution to the wave equation. The constant A in Eq. 
(4. 19) combines the constants Ay and B r of the preceding equations. Note that we have 
replaced the approximation symbol with the equality sign in Eq. (4.19). This is only a 
matter of notation, and we still recognize that the solution is only an approximation. 

Since the linear wave equation is in terms of partial differentials with respect to y and r, 
if we multiply Eq. (4.19) by a function of k r , then the result is also a solution to the wave 
equation since k r is independent of y and r [Ref. 9:p. 125]. Similarly, any summation of 
such solutions will also be a solution to the wave equation. Combining these concepts into 
a generalized integral form results in 



where A(k r ) is the term that incorporates the constant A of Eq. (4.19) and the arbitrary 
function of k r as discussed. Equation (4.21) is therefore our solution of the linear wave 
equation based on the WKB solution for an inhomogeneous medium. If Eq. (4.21) is 
correct, then for the case of a homogeneous medium it should collapse into the free-space 
Green’s function for an isospeed medium. As well as validating the form of the equation, 
it should also provide us with an expression for the term A(k r ). 




(4.21) 
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E. THE FREE-SPACE GREEN'S FUNCTION 

The free-space problem for an isospeed medium is one in which there are no 
boundaries to consider and the medium is considered to be homogeneous. Since one of 
the assumptions made in the development of Eq. (4.21) was that there was no surface or 
bottom boundary interaction, all that we have to do to find a solution for the free-space, 
isospeed problem is to apply the conditions for a homogeneous medium to Eq. (4.21). 
For the homogeneous medium, k Y (y) = k Y (y 0 ) is a constant. For this case, there is a 
closed form solution of the integral appearing in the exponent of Eq. (4.21), which results 
in a solution of the form 



It should be noted that for the case of a homogeneous medium, the simplifying assumption 
given by Eq. (4.17) is unnecessary, and Eq. (4.22) is an exact solution. 

We will now specifically solve the time-independent Helmholtz equation of the 
form of Eq. (4.3) for the free-space, isospeed problem and compare the result with Eq. 

(4.22). In the free space problem a wave will radiate in an outwards direction, 
experiencing spherical spreading, and be free of boundary interactions. The homogeneous 
medium condition means that k(r) = k is a constant. Since the wave will spread 
spherically, it makes sense to use the spherical coordinate system to describe the solution. 
It can be shown that 




(4.22) 




(4.23) 
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is a solution to Eq. (4.3), under the conditions described, everywhere except at r = 0. 
Note that in Eq. (4.23) r is a scalar quantity rather than a vector. Equation (4.23), 
however, does satisfy 



for all r. This is consistent with there having to be a source at r = 0 to generate the wave. 
The choice of a point source as the source is a convenient one since in Chapter 2 we 
assumed that our transmit aperture was an array of point sources. If we think of this 
forcing function or source term as having unit amplitude, we set the constant in Eq. (4.23) 
equal to (-1/471). Applying this amplitude term and converting to cylindrical coordinates 
we can show that 



V 2 tp(r) + k 2 <p(r) = -4jcA5(r) 



(4.24) 




(4.25) 



is a solution to 



V 2 cp(r) + k 2 <p(r) = 8(y-y 0 ) 



(4.26) 



for all r where 




(4.27) 
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and r is the polar radial distance. Equations (4.23) and (4.25) are the free-space Green's 
functions for an isospeed medium in spherical and cylindrical coordinates, respectively. 
They represent the solution to the time-independent Helmholtz equation for a homogeneous 
medium with no boundary interactions, and matching the boundary condition for a point 
source. 

Officer [Ref. 9:p. 126] shows that the free-space Green's function in cylindrical 
coordinates can be expressed as an integral of the zero-order Bessel function of the first 
kind, having the form 



^jkR 

R 



/; 



k ’ j. n . rte jkY(5 '° )(y ‘ y ° , dk 

jky(y 0 ) 0( r ) r - 



(4.28) 



If we now combine Eqs. (4.28) and (4.25) we obtain 

<P (r) =J^ 4n ky'typ) J ° <k r r)e '' llY(y ° >(y y ° >dk r (4.29) 

which is the solution found by solving the linear wave equation specifically for the free- 
space, isospeed problem. 

F. DETERMINATION OF THE OCEAN MEDIUM TRANSFER FUNCTION 

Returning to our solution for the free-space problem based on the WKB approximation 

given by Eq. (4.22), we see that it is in a similar form to Eq. (4.29) but is in terms of a 

Hankel function rather than a Bessel function. Based on observations made earlier in this 

chapter, however, we would expect to be able to express Eq. (4.22) in terms of a Bessel 

2 1 

function. This is in fact the case, and using the relationships Ho(-k r r) = -Ho(k r r) and 
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J o (k r r) = j [ Ho(k r r) + H(j(k r r) ] from Brekhovskikh [Ref. 8:pp. 76-77] it can be 



shown that 



<p(r) = 



0 



2A(k r ) Jp(k r r) ^.j k Y (y 0 ) (y-y c ) 

V|k Y (y 0 )| 



(4.30) 



Our two solutions, given by Eqs. (4.29) and (4.30), are now in a comparable form and 
we can determine the value of the constant A(k r ) to be 



Note the fact that the two terms involving the wave number in the vertical or depth 
direction, k Y (y 0 ), in Eq. (4.31) are not partially cancelled. This is due to the fact that 
their value may be either purely real (propagating wave) or purely imaginary (evanescent 
wave) caused by the infinite range of k r . As a consequence, Eq. (4.31) is the most 
concise way of expressing this quantity. 

The fact that we were able to get our solution from the WKB approximation, Eq. 
(4.30), into the same form as the Green’s function, Eq. (4.29), for the free-space,isospeed 
problem suggests that our form of the WKB solution is correct. Combining Eqs. (4.31), 
(4.21), and changing the variable of integration from k r to f r it can be shown that 




(4.31) 
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which is our solution based on the WKB approximation for a sound-speed profile that is a 
function of depth only. This is one of the fundamental results of this chapter. Now that 
this quantity is known, we can calculate the ocean medium transfer function based on the 
WKB approximation. 

The simplest way to calculate the ocean medium transfer function is to compare Eq. 
(4.32) with the overall system complex frequency response of Eq. (3.14) calculated for a 
single point source of unit amplitude. Under these circumstances MT, NT, and c mn (f) 
equal unity and Eq. (3.14) reduces to 

H(f,r) = 2it f H M (f,f r ,y o ;y)J 0 (2itf r r m )e' i2 ’ tfY(y « K) '' y ” ) f r <lf r . (4.33) 
J o 



The reason that we are interested in Eq. (4.33) is that these are the very same conditions 
that led to the solution given by Eq. (4.32). Consequently, Eqs. (4.33) and (4.32) should 
be the same. Any differences between these two equations must therefore be accounted 
for by the ocean medium transfer function as this is the only unknown quantity in the 
equations. It can therefore be shown that 



H M (f,f r ,y 0 ;y)= 

2k Y (y 0 Mk Y (y)| 



^yo 



where k Y (y) is as defined by Eqs. (4.20a) and (4.20b). Equation (4.34), then, is our 
ocean medium transfer function based on the WKB approximation and is the fundamental 
result of this chapter. 
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V. THE SIGNAL GENERATOR 



A. THE COMPLEX ENVELOPE 

Consider a transmitted electrical signal of the following form : 

x(t) = a(t)cos( 2jtf c t + 0(t) ) (5.1) 

This is an amplitude and angle modulated carrier where we will assume that a(t) and @(t) 
are arbitrary, real baseband amplitude and angle modulating signals, respectively. This 
implies that both a(t) and 0(t) are bandlimited to the region -W < f < W. These signal 
components, a(t) and 0(t), are commonly known as the natural envelope and phase of the 
signal x(t), respectively. As a consequence, x(t) is a real bandpass or bandlimited 
signal. This means that the amplitude spectrum of the signal is zero outside a frequency 
band of width 2W Hz centered about +/- f c Hz, where f c > W Hz and f c is referred to as 
the center or carrier frequency. 

Using Eq. (5.1) for the model of x(t), we will make use of the concept of representing 
a real bandpass signal in terms of its complex envelope [Ref. 3:pp.l76-188] [Ref. 10:pp. 
74-90]. The complex envelope of x(t), denoted by x(t), is defined as 

x(t) = [ x(t) + jx(t) ] e' j27cfct (5.2) 

where x(t) is the Hilbert transform of x(t). If x(t) is real, then its Fourier transform 
exhibits Hermitian symmetry. Consequently, the frequency spectrum along the negative 
frequency axis is completely determined by the spectrum for positive frequencies (and visa- 
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versa). Thus, if we remove the frequency spectrum along the negative (or positive) 
frequency axis, we will not loose any of the information contained in the signal. 

The Hilbert transform of x(t) is defined by 



which is the convolution of x(t) with the time function l/7tt. This has the effect of 
producing a -/+ 90° phase shift for all positive/negative frequencies of the input signal. 
The Hilbert transform applies to any signal that is Fourier transformable, and by looking at 
Eq. (5.3) in the frequency domain, we obtain the analogous definition 



where sgn(f) is the signum function and IFT{*} is the inverse Fourier transform. The 
effect of adding a real signal to its Hilbert transform is to remove the negative frequency 
components, producing a one-sided frequency spectrum which is known as the 
pre-envelope x+ (t). [Ref. 10:pp. 74-75] [Ref.ll:pp. 67-72] 

This one-sided spectrum is then translated down in frequency so as to be centered at 0 
Hz to obtain the complex envelope. Thus, the complex envelope of a real bandpass signal 
is typically a complex signal with a lowpass or baseband frequency spectrum that has the 
same information content as the original bandpass signal. 

Applying the definition of the Hilbert transform to the specific case of the real bandpass 
signal x(t) in the form of Eq. (5.1) results in 




(5-3) 



x(t) = IFT { -jsgn(f) X(f) } 



(5.4) 
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x(t) = X c (t)sin(27cf c t) + x s (t)cos(2jtf c t) 



(5.5a) 



where 



x c (t) = a(t)cos0(t). 


(5.5b) 


x s (t) = a(t)sin0(t). 


(5.5c) 



and x c (t) and x s (t) are known as the in-phase and quadrature components of the bandpass 
signal, respectively. It should be reiterated here that one of the underlying assumptions in 
this treatment of Eq. (5.1) is that x(t) is a bandpass signal. This in turn requires x c (t) and 
x s (t) to be bandlimited to the interval -W < f < W where f c > W. Equation (5.5a) is 
sometimes referred to as the canonical form of x(t). [Ref. 3:p. 183] [Ref. 10:p. 85] 

By combining Eqs. (5.1), (5.2) and (5.5), it can be shown that 

x(t) = a(t)e j0(t) (5.6) 

is the complex envelope of the assumed form of the transmitted electrical signal given by 
Eq. (5.1). Given our initial assumptions, this means that the amplitude spectrum of the 
complex envelope described by Eq. (5.6) will be baseband centered around 0 Hz with a 
bandwidth of W Hz. 

The form of Eq. (5.6) demonstrates the advantages of complex envelope notation. 
First, it allows for the simple representation of amplitude and angle modulation. Second, 
for our assumed form of x(t) given by Eq. (5.1), which was a modulated carrier, the 
complex envelope is independent of the carrier frequency. The significance of this 



48 



independence from the carrier becomes apparent when dealing with the discrete 
representation of x(t), required for signal processing by digital computers. Since the 
information content of the signal x(t) is completely represented by its complex envelope 
x(t), we can sample the baseband complex envelope rather than the original bandpass 
signal. As the maximum frequency component present in x(t) is W Hz rather than 
(f c + W) Hz as in x(t), a lower sampling rate is required to satisfy Nyquist's criterion. 
Consequently, less time samples or data points will be required to represent the signal 
resulting in a reduction of computation time. 

It can be shown that both the original bandpass signal x(t) and its Fourier transform 
X(f) can be expressed in terms of the complex envelope using 



From the previous equations it can be seen that it is a straight forward matter to move 
between the complex envelope and common signal representations of x(t) in both the time 
and frequency domains. Thus the idea of working with the complex envelope is not at the 
expense of any substantial computational penalty in going from one representation of the 
signal to the other. [Ref. 10:pp. 85-86] 



x(t) = Re[x(t)e^ 2,rfct ] 



(5.7) 



and 




(5.8) 
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B. THE CW AND LFM PULSE 



For the purposes of this thesis, two classic waveforms used in underwater acoustics 
will be modelled. They are the CW (Continuous Wave) and LFM (Linear Frequency 
Modulated) pulses, both of which are real bandpass signals that can easily be expressed in 
the general form of Eq. (5.1). 

A CW waveform is simply one in which there is no angle modulation. Thus, the 
expressions for the bandpass CW signal and its complex envelope can easily be found by 
setting 0(t) equal to zero in Eqs. (5.1) and (5.6) yielding 



Frequency Modulation is one of the most commonly used types of angle modulation 
and is where the instantaneous frequency, fi (t), is varied linearly with some arbitrary 
baseband signal. The instantaneous frequency is related to the time derivative of the 
instantaneous phase and, for the real bandpass signal of the form of Eq. (5.1), it can be 
shown that 



x cw(0 = a(t)cos(2itf c t), 



(5.9) 



and 



x cw (0 — a(t). 



(5.10) 



f i (t) = f c + _U 



J_ d(effl) 

2 k dt 



(5.11) 



where f c is the carrier frequency. [Ref. 10:pp. 180-183] 
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If the instantaneous frequency is a linear function of time rather than some arbitrary 
baseband signal, then Eq. (5.1) is referred to as a LFM waveform. Let the angle 



modulating signal be of the form 



0(t) = D p t' 



2 



(5.12) 



where Dp is referred to as the phase-deviation constant or phase sensitivity. Substituting 
Eq. (5.12) into Eq. (5.11) results in the following expression for the instantaneous 
frequency : 



The value of the instantaneous frequency given in Eq. (5.13) is a linear function of time as 
desired, composed of a constant carrier frequency and the term (Dpt / 7t) which is known as 
the frequency deviation. The value (Dp / k ) in this instance is sometimes referred to as the 
frequency sensitivity of the modulator. [Ref. 3:pp. 193-200] [Ref. 12:pp. 267-271] 

Upon substituting Eq. (5.12) into Eqs. (5.1) and (5.6), we obtain the following 
expressions for the LFM waveform and its complex envelope : 




(5.13) 



x lfm (0 — a(t)cos(2jtf c t + D p t ) 



(5.14) 



and 




(5.15) 
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On comparing the above expressions with Eqs. (5.9) and (5.10), it can be seen that the 
expressions for the CW waveform can be easily obtained from those for the LFM 
waveform by setting the phase deviation constant equal to zero. This has important 
consequences for the computer code that generates these waveforms. It means that we 
only require code for the LFM waveform since we have seen that CW is merely a special 
case of LFM. 

For the purposes of this thesis we will be concerned primarily with the rectangular 
amplitude modulating function given by 



a(t) = { 



A , |t|ST p 

0 , |t| > Tp 



(5.16) 



where A is a constant and Tp is the pulse length (in seconds). It should be noted from Eq. 
(5.16) that the pulse is assumed to be centered at t = 0, and that the modulating function is 
in essence acting like a weighted on/off switch, either scaling the waveform by a constant 
or setting it equal to zero. The reason that we are interested in this particular amplitude 
modulating function is that it is because it corresponds to that transmitted by a "typical" 
SONAR system. Other amplitude modulating functions that will be considered are the 
Hanning and Hamming weighting functions, which will be discussed in the final section of 
this chapter. 

For a LFM pulse with a pulse length of Tp (in seconds), the quantity (D p T p / 7t) is 
referred to as the swept bandwidth (SBW) of the signal x(t). This value represents the 
frequency change (in Hz) exhibited by the signal over the period of a single pulse and 
assumes that only one sweep occurs during the period of the pulse. We are further 
assuming that this swept bandwidth is centered around the carrier frequency. 
[Ref. 3:pp. 197] 
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C. FOURIER SERIES REPRESENTATION 



As we are modelling the ocean medium in terms of transfer functions, it is easier and 
more time efficient to work in the frequency domain as opposed to the time domain. 
Accordingly, we require x(t) to be in a form for which the transform X(f) can be easily 
computed. Except for a few simple cases, there is no closed form expression for X(f) . In 
fact, X(f) for the rectangular envelope LFM pulse is similar in form to a Fresnel diffraction 
integral. 

Assume that we can expand our complex envelope into the truncated complex Fourier 
series 



x(t) = 



£ c q e W »' 
q=-K 




(5.17) 



where 



c 



q 



To 



,T 0 /2 
•'-T 0 / 2 



x(t)e 



-j2irqf o t 



dt 



(5.18) 



is the complex Fourier series coefficient at harmonic q, f 0 = (1/T 0 ), and T 0 is sometimes 
referred to as the data record length and is the time interval over which we are attempting to 
approximate the signal. Since this series is expanded about a period of T 0 , the 
fundamental frequency of the Fourier series expansion is f 0 . Accordingly, the coefficients 
represent the magnitude and phase of the frequency components at multiples of f 0 Hz. 
This can be clearly demonstrated by taking the Fourier transform of Eq. (5.17) which 
yields the result: 
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(5.19) 



K 

X(f)= Z c q 8(f-qf 0 ). 
q=-K 

The value K determines the number of terms contained in the truncated Fourier series and 
can be interpreted as the number of harmonics required to represent the complex envelope 
x(t). 

The effectiveness of the series representation given by Eq. (5.17) lies in the ability of a 
small number of terms to provide an adequate approximation of the original signal. The 
complex exponentials 



used in this series expansion form an orthogonal basis set of functions [Ref. 15:pp. 214- 
215]. Consequently, the approximation from Eq. (5.17) using the coefficients calculated 
from Eq. (5.18) results in the minimum mean squared error estimate of x(t) for the number 
of terms used [Ref. 1 1 :pp. 32-36]. The choice of complex exponentials as the basis set of 
functions produces a series expansion that relates the time and frequency domains. We 
have to consider a truncated series because we can only process a finite amount of data, and 
the less terms we need to represent x(t) the less computer processing will be required. 

The preceding discussion of the Fourier series assumed that x(t) was a continuous time 
signal but the concepts and results can equally be applied to a discrete time version. We 
will now discretize x(t) and introduce the Discrete Fourier Transform ( DFT{-}). The 
forward and inverse DFT can be defined as 
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(5.21) 



DFT(x(l))= £ x(De' i2,Rll/L 
1 = -L' 

and 

IDFT(X(q))=i £ X(q)e j2 " ql/L , (5.22) 

q = -K 

respectively, where IDFT{ • } is the Inverse Discrete Fourier Transform and L = (2L* + 1) 
is the total number of time samples taken of x(t) during the time period T 0 . To digitize 
x(t), we will sample it at time intervals of (T 0 /L) which results in x(l) = x(lT 0 /L). 
Substituting this digitized version of x(t) into Eqs. (5.17) through (5.19) and comparing 
with Eqs. (5.21) and (5.22), it can be shown that 

c,=iX(q)=iDFTWl)) (5.23) 



and 



x(l) =LxIDFT{c q ). (5.24) 

Equations (5.23) and (5.24) in conjunction with Eqs (5.7) and (5.8) are the primary 
equations that we will use in moving between the time and frequency domains of both the 
bandpass signal and its complex envelope. 

It should be noted that this concept of representing a signal by a truncated Fourier series 
is completely general and can be applied to any signal. No assumptions regarding the 
form of the complex envelope, x(t), were made in the development of this section. 
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Consequently, if the complex envelope is of the form given by Eq. (5.6), no assumptions 
have been made regarding the amplitude and angle modulation terms of Eq. (5.1). 

D. IMPLEMENTATION 

Let us now consider the specific case of rectangular amplitude modulated CW and LFM 
pulses. We require a method of determining the value of K, so that the Fourier series 
representation contains an adequate number of harmonics to accurately represent x(t), and 
a value of L such that the sampling interval is sufficient to satisfy Nyquist’s criterion or 
some other specified sampling rate. 

The complex envelope of a CW pulse can be thought of in the time domain as a 
rectangular window convolved with a sequence of delta functions at intervals of T c 
seconds. It can easily be seen that as a consequence, the resulting DFT will be a discrete 
representation of the sine waveform. To represent this signal we will take all the 
frequency components (separated by 1/T 0 Hz) occurring within three sine function zero 
crossings (each separated by 1/Tp Hz) from the origin. All the frequency components 
extending past this point will be less than 10% of the value of the frequency component at 
the origin, so we will consider them as being insignificant and ignore them, which yields 
the relationship 



K 



cw 




(5.25) 



As stated earlier, there is no closed form expression for the LFM pulse in the frequency 
domain. However, for a rectangular envelope LFM pulse Papoulis [Ref. 1 2:pp. 267-27 1 ] 
shows that if the following condition is satisfied 



56 



(2f c -SBW)»Y^ , 



(5.26) 



then the resulting transform can be approximated as having a constant magnitude within a 
frequency band equal to the swept bandwidth, centered about the carrier frequency f c . 
Outside of this region the transform decays to zero. We need to approximate this rate of 
decay so as to decide how many frequency components we need to accurately represent the 
original bandpass signal x(t). 

The corresponding complex envelope of the LFM pulse will therefore be approximately 
constant within the region -SBW/2 < f < SBW/2 and decaying to zero outside of this 
region. Since the frequency components are spaced 1/T 0 Hz apart, the value of K required 
to represent this constant region of the spectrum is given by 



K 



SBW T 0 
2 



(5.27) 



Now consider the CW pulse from the point of view as being a special case of the LFM 
pulse where the swept bandwidth equals zero. We would therefore expect that in the limit 
as the swept bandwidth goes to zero that Kif m = Kc W . We will therefore approximate 
this decaying region of the spectrum as being the width of three zero crossings of the sine 
function resulting from just the rectangular amplitude function. It is surmised that 
frequency components outside this region will be less than 10% of the frequency 
component at the origin and can therefore be disregarded. Therefore, adding the K values 
from Eqs. (5.25) and (5.27) yields 



SBW T 0 
2 



3T P 

T n 



(5.28) 
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which is valid for both CW and LFM pulses. It should be noted from Eq. (5.17) that K is 
an integer value, so in our computer implementation of Eq. (5.28), K is rounded up to the 
next highest integer so as not to increase the truncation error. 

Like the value for K, L is also an integer, but with the additional constraint that it is an 
odd number. This is in order to satisfy the assumed form of Eq. (5.21) and to make the 
simulation output conform to that required by the adaptive beamforming algorithm that will 
be used to validate some of our results. The number of samples required can easily be 
computed from the calculated values of K using 

L = SR x K (5.29) 

where SR is the user defined value of the sampling rate and K is the value determined 
from Eq. (5.28). In order to satisfy Nyquist's criterion, the sampling rate must be greater 
than two. This means that the sampling frequency must be greater than twice the highest 
frequency component present in the truncated series representation of the transmitted 
electrical signal. In the computer implementation of Eq. (5.29), L is rounded to the next 
highest odd integer value in order to satisfy all of the previously mentioned constraints. 

In the computer code, the amplitude and angle modulated signal of Eq. (5.1) is defined 
by the following variables whose values are supplied by the user: 

• AMP : the magnitude of the rectangular amplitude modulating signal. 

• TP : the pulse length of the transmitted pulse (sec). 

• PRF : the pulse repetition frequency (Hz). This value is the same as f Q and is 

equal to the inverse of the data record length T 0 (sec). 

• FC : the carrier frequency (Hz). 

• SBW : the swept bandwidth (Hz). To describe a CW pulse, set this 

input parameter to zero. 
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• CHIRP : indicates whether the LFM sweep is increasing in frequency (up) or 

decreasing (down). 

• SR : sampling rate. This value is the ratio of the sampling frequency to the 

maximum frequency component present in the signal and must be greater 
than two in order to satisfy Nyquist's criterion. 

• AMOD : the amplitude modulation function, selected from a menu of choices. The 

rectangular amplitude modulation function is denoted by AMOD = 0. 

E. THE RECTANGULAR AMPLITUDE MODULATED PULSE 

To validate the theory of the preceding section, examples of rectangular amplitude 
modulated CW and LFM pulses were reconstructed from the calculated Fourier series 
coefficients of its complex envelope, denoted xc e (n), and compared with the actual 
sampled bandpass signal denoted x(n). In the following plots, x(n) is plotted as a dashed 
curve while xc e (n) is plotted as a solid curve. For plotting purposes, a third-degree 
polynomial curve fitting routine is used to interpolate between the data points. 

1. CVV Pulse with 50% Duty Cycle 

Duty cycle is defined as the ratio Tp/T 0 , which for this case means that the 
amplitude weighting function is nonzero for 50% of the time. While this is not a 
particularly practical transmit signal for underwater acoustic applications, it will 
demonstrate that a real bandpass signal can be reconstructed from its complex envelope. 
Combining this high value of duty cycle with a suitably low carrier frequency, the 
comparative plots of x ce (n) and x(n) will be uncluttered and will help illustrate some 
important points. 

a . Pulse Reconstruction 

A CW pulse with the following characteristics is shown in Fig. 5.1a 
• AMP = 1, TP = 0.0125 sec, PRF = 40 Hz, FC = 400 Hz, SR = 4.0 
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It can be seen from Fig. 5.1a that the two curves, x(n) and x ce (n), agree 
fairly closely. This indicates that we are indeed able to reconstruct a real bandpass signal 
from its complex envelope. However, there are some important effects to be observed. 
Note the manner in which Xc e (n) overshoots x(n) when the pulse is turned on and off. It 
can also be seen that Xc e (n) tends to oscillate around the true value of x(n). 

A truncated Fourier series, by its very nature, will exhibit high frequency 
oscillations about the true value of the function. While these oscillations usually have 
sufficiently small amplitudes to be of no major consequence, this is not the case at a 
discontinuity. Where x(n) is turned on/off corresponds to just such a discontinuity in the 
rectangular amplitude weighting function. As more terms are added to the truncated series 
the overshoot does not tend to disappear. Instead, it becomes narrower due to the higher 
frequency terms present and the first overshoot approaches a value of 0.08949 of the 
discontinuity. This behavior is known as "Gibbs Phenomenon" and is due to the truncated 
Fourier series representation of the complex envelope, rather than the fact that we have 
used the complex envelope to reconstruct the original bandpass signal. 
[Ref. 13:pp. 107-112] [Ref. 14:pp. 179-185] 
b. Lanczos Smoothing 

It was suggested by Cornelius Lanczos [Ref. 15:pp. 219-221] that this 
phenomenon could be reduced if the Fourier series could be made to converge more 
quickly. He showed that if the series coefficients were modified by a correction term he 
called the sigma factors, given by 



o(K,q) = 



sin(qrc/K) 
(qrc/K) . 



(5.30) 
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where K is the order of the last term in the series, then the tendency of the high order terms 
to make a series divergent could be counteracted. As can be seen from Eq. (5.30) this was 
achieved by increasing the amount of attenuation with increasing order of the Fourier 
coefficients. Hamming [Ref. 13:pp. 107-112] comments that K in Eq. (5.30) is usually 
replaced with (K+l), and this is the convention we will adopt in our application of the 
sigma factors for this thesis. 

Lanczos also showed that the application of the sigma factors was 
analogous to averaging over an incremental time period of ±T 0 / 2K. This means that each 
original point in the truncated Fourier series, x(t), is replaced by the arithmetic mean around 
that point, between the limits of the incremental time period. Since the ripple in the 
truncated series representation has the period of either the first term neglected or the last 
term kept, the main effects of this ripple can be removed by integrating or averaging over 
this period. The averaging interval usually chosen for this smoothing is T 0 /(K+l), which 
corresponds to replacing K with (K+l) in Eq. (5.30). If K is large, then the region of 
local averaging is very small and has very little effect, except in the neighborhood of the 
discontinuity. This is a desirable result since the discontinuity is the region where the 
smoothing is most needed. 

When applied to the truncated Fourier series, this is known as "Lanczos 
Smoothing". Applying this smoothing to the truncated series of Eq. (5.17) yields the 
following: 



*ls(t) : 



: £ G(K,q) Cq e J 

q=-K 



j2rajf 0 t 



^ T 0 
"T 



(5-31) 
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As can be seen from Eq. (5.31), the smoothed series is nothing more than the original 
Fourier series with its coefficients multiplied by the corresponding sigma factors. This 
result is what we would expect owing to the linearity of the Fourier series. The effect of 
time averaging is the same as convolving the time sequence with a weighted rectangular 
window. This in turn is analogous to multiplying the corresponding Fourier series 
coefficients of the time sequence with those of the window. Thus, the sigma factors can 
be thought of as the Fourier series coefficients of this weighted rectangular window, 
c. The Smoothed Pulse 

The effect of the application of these sigma factors can be seen in Fig. 5.1b, 
where the Lanczos Smoothing has been applied to the Fourier Series coefficients of our 
CW pulse. As expected, the oscillations and overshoot at the discontinuity have been 
reduced but at the expense of slightly broadening the pulse and reducing the rate of change 
at the leading and trailing edges of the pulse. It should be noticed that the oscillations have 
been reduced for both portions of the duty cycle, i.e., not just during the pulse. 

The first effect is to be expected since we are effectively convolving our 
time sequence with a windowing function. Whenever we convolve two sequences the 
result is a sequence larger than, but related to, the size of the original functions. Since the 
windowing function is significantly smaller than the pulse duration, we expect the effects 
of this spreading to be minimal. On comparing Figs. 5.1a and 5.1b we can see that this is 
certainly the case for this example. 

The second effect is due to the fact that the smoothing technique attenuates 
the high order Fourier coefficients. Attenuating high order coefficients is analogous to 
attenuating the high frequency components of the series representation of the bandpass 
signal. Thus, we would expect the smoothed series representation of the signal, Xc e (n), to 
be slower than the original Fourier series at reacting to sudden change, i.e., a discontinuity. 
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CASE: FIG 5.IR WRVEFORM:CW PULSE ( NFREO: 13 ) 
fi: 1.0 IP: 12.5 MSEC PRF: <0.0 HZ FC: <00.0 HZ 



Fig. 5.1a C\V pulse, 50% duty cycle, no smoothing. 




CRSE--FIG 5. IB WnVEFORM : CW TUL5E ( NFREO: 13 ) 

fl: 1.0 TP: 12.5 MSEC FRF: <0.0 HZ FC: <00.0 HZ 



Fig. 5.1b CW pulse, 50% duty cycle, Lanczos smoothing. 
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Regarding the application of the sigma factors to the specific problem of Gibbs 
phenomenon, Lanczos [Ref. 15:pp. 225-226] notes that unlike the often-used Fej6r 
arithmetic mean method, the sigma factors do not completely eliminate the oscillations due 
to the discontinuity. However, the oscillations are significantly reduced and the sigma 
factors do not cause nearly as drastic a reduction in the rate of change at the leading and 
trailing edges of the pulse. It should be noted that this smoothing invalidates the 
minimum-mean-square-error property of the Fourier series. The original Fourier series 
coefficients are the only ones that satisfy this property and now they have been changed. 

2. LFM Pulse with 50% Duty Cycle 

A LFM pulse with the following characteristics is shown in Fig. 5.2a 



• AMP = 1, TP = 0.0125 sec, PRF = 40 Hz, FC = 400 Hz, SRATE = 4.0 
SBW = 60Hz, "up" chirp 



The purpose of Figs 5.2a through 5.2c is to demonstrate that our analysis of the LFM pulse 
was correct. As can be seen from Fig. 5.2a, we are able to reconstruct the LFM 
waveform from the Fourier series representation of its complex envelope, as was the case 
for the CW pulse. The LFM pulse in this case is said to have "up chirp" because its phase 
deviation constant Dp is positive, that is, the frequency is increasing during the period of 
the pulse. 

As for the CW pulse, the reconstructed signal x ce (n) oscillates about its true 
value x(n) and exhibits the characteristic overshoots at the leading and trailing edges of the 
pulse. The application of the sigma factors in Fig. 5.2b shows that the Lanczos smoothing 
is as effective for the LFM case as it was for the CW pulse and exhibits the same results. 
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TRANSMITTED ELECTRICAL SIGNAL X 




Fig. 5.2a LFM pulse, 50% duty cycle, "up chirp", no smoothing. 




A: 1.0 TP: 12.5 MSEC PRF : <0.0 HZ PC: <00.0 HZ SWPIBW: 60.0 HZ CHIRP: U 

Fig. 5.2b LFM pulse, 50% duty cycle, "up chirp", Lanczos smoothing. 
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It sh6uld be noted however, that since the swept bandwidth is low and the pulse 
repetition frequency is fairly high, the value of K given by Eq.(5.27) is less than one. 
This effectively means that the signal is approaching the limit of becoming a CW pulse. 
Inspection of Eq. (5.28) shows that the dominating effect on the number of significant 
frequency components is due to the rectangular pulsed nature of the waveform rather than 
the frequency modulation. ' Further investigation needs to be carried out on whether or not 
Eq. (5.27) is a valid approximation. This will be done by the study of a LFM pulse with a 
smaller duty cycle to follow. 

The LFM pulse in Fig. 5.2c differs from the previous LFM pulses in that it is a 
"down chirp" pulse, that is its phase deviation constant is negative and the frequency is 
decreasing during the period of the pulse. The purpose of this Figure is to show that the 
computer code can model the "down chirp" pulse as well as the "up chirp" case. 




TIME T(MSEC) 

CRSEs FIG 5.2C WAVEFORM : LFM PULSE ( NFREQ: 15 ) 

fl: 1.0 TPs 12.5 MSEC PRFs 40.0 HZ FC: tOO.O HZ SWPTBW: 60.0 HZ CHIRP: D 



Fig. 5.2c LFM pulse, 50% duty cycle, "down chirp", Lanczos smoothing. 
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3. CW PULSE WITH 8% DUTY CYCLE 



A pulse waveform with a small duty cycle is of more practical use for 
underwater acoustics applications, so we need to check that our computer code will handle 
more realistic signals. A CW pulse with the following characteristics is shown in Figs. 
5.3a and 5.3b: 



• AMP = 1, TP = 0.04 sec, PRF = 2.0 Hz, FC = 400 Hz, SR = 4.0 



As can be seen from the above parameters, the ratio T p /T 0 yields the value 0.08 which 
corresponds to a duty cycle of 8%. In Figs. 5.3a and 5.3b, only the curve reconstructed 
from its Fourier series coefficients, xce(n), is shown so that the plots do not appear too 
cluttered. 

From Fig. 5.3a we can observe that as expected, the unmodified Fourier 
coefficients produce a waveform that oscillates about the true value of x(n). The ringing, 
or oscillation, at the leading and trailing edges of the pulse is quite pronounced and give the 
appearance of the rectangular pulse shape spreading. 

The application of Lanczos smoothing to this case, shown in Fig. 5.3b, 
demonstrates that the oscillations about the true value have been reduced but not completely 
removed. This can be seen to be at the expense of a slight spreading of the pulse and a 
decrease in the rate of change at the leading and trailing edges of the pulse. As expected, 
these are the same effects as were observed for the CW pulse with the 50% duty cycle 
except that the extent of the ringing present in this case is more dramatic. 
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R: 1.0 TPs 10.0 MSEC PRF: 2.0 HZ FC: •100.0 HZ 

• Fig. 5.3a CW pulse, 8% duty cycle, no smoothing. 




fl: 1.0 TP: 40.0 MSEC PRF: 2.0 HZ FC: 400.0 HZ 

Fig. 5.3b CW pulse, 8% duty cycle, Lanczos smoothing. 



68 



Why are we interested in this amount of ringing in the pulse representation? 
One of the outputs from this computer simulation will be the time-domain output electrical 
signal at each element in the receive array. We wish to use this information to illustrate the 
pulse distortion due to dispersion in the ocean waveguide. Consequendy, we wish to 
"clean up" the transmitted pulse as much as possible so that we do not confuse distortion 
due to the Fourier series representation of the transmitted signal with distortion due to 
dispersion in the ocean waveguide. 

By comparing Figs. 5.3a and 5.3b, it can be seen that the application of 
Lanczos smoothing does produce a much "cleaner" pulse shape. It is our contention that it 
will be easier to observe the effects of dispersion on this smoothed pulse that on the 
unmodified pulse. 

4. LFM Pulse With 8% Duty Cycle 

A LFM pulse with the following characteristics is shown in Figs. 5.4a and 

5.4b: 



AMP = 1, TP = 0.04 sec, PRF = 2.0 Hz, FC = 400 Hz, SRATE = 4.0 
SBW = 120 Hz, "up chirp" 



As for the CW case, we have only presented the signal reconstructed from its Fourier series 
coefficients, xce(n), so that the plots do not become too cluttered. 

It can be observed that all the comments made regarding the CW pulse with the 
8% duty cycle are applicable here. It can be seen from comparing Figs. 5.3a and 5.4a 
however, that the oscillations about the true value and the amount of ringing appears to be 
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CRSE: FIG 5.40 WnVEf'ORfl ! LF'M PULSE ( NFREQ: 137 ) 

R : 1.0 TP: 40.0 MSEC FRF: 2.0 HZ FC: 400.0 HZ SWFTBW: 120.0 HZ CHIRP:? 



Fig. 5.4a LFM pulse, 8% duty cycle, "up chirp", no smoothing. 




R: 1.0 TP: 40.0 MSEC FRF: 2.0 HZ FC: 400.0 HZ SWPTBW: 120.0 HZ CHIRP: 7 

Fig. 5.4b LFM pulse, 8% duty cycle, "up chirp", Lanczos smoothing. 
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less for the LFM pulse than it does for the CW pulse. This is probably due to the number 
of frequency components used to represent the signal. 



It would appear that the value of K given by Eq. (5.28) is a more conservative 
estimate of the number of frequency components required to represent the LFM pulse than 
the estimation for the CW pulse given by Eq. (5.25). As the number of frequency 
components used is increased, the amount of ringing at the leading and trailing edges 
should be reduced. Since the amount of ringing for the LFM waveform is less than for the 
CW case, we assume that we erred on the side of caution and included more frequency 
components that absolutely necessary to represent a recognizable LFM pulse shape. This 

is consistent with the logic we used in proposing the form of Eq. (5.28) 

This would also explain why the application of the Lanczos smoothing in Fig. 
5.4b appears to have distorted the LFM pulse shape slightly less than it did for the CW 
pulse in Fig. 5.3b. If we are at the minimum number of frequency components required to 
represent the signal, then the attenuation of the high order terms will have a more dramatic 
effect than if we have a slight excess in the number of frequency components required. 
This effect can be shown quite dramatically in Figs. 5.4c and 5.4d. 

In Figs. 5.4c and 5.4d, we have reconstructed the LFM pulse using the number 
of frequency components given by Eq. (5.27), which assumes that only the constant 
portion of the spectrum of the LFM pulse is significant. As can be seen from Fig. 5.4c, 
the reconstructed pulse xce(n) is far from being a clean looking pulse shape. It oscillates 
erratically about the true value and the amount of ringing exhibited is quite extensive. 
When we apply the Lanczos smoothing to this case the waveform appears distorted as 
expected. In fact, it is distorted to the extent that we are not able to tell by looking at it that 
it has a rectangular pulse shape. 
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R: 1.0 TP s 40.0 MSEC PRF: 2.0 HZ PC: 400.0 HZ SWPTBW: 120.0 HZ CHIRP: U 

Fig. 5.4c LFM pulse, 8% duty cycle, "up chirp", no smoothing, 
insufficient frequency components. 




CASE: FIG 5.40 WRVEFORMrLFM PULSE ( NFREQ: 63 1 

R: 1.0 TP: 40.0 MSEC PRF: 2.0 HZ FC: 400.0 HZ 5WPTBW: 120.0 HZ CHIRP: U 



Fig. 5.4d LFM pulse, 8% duty cycle, "up chirp", Lanczos smoothing, 
insufficient frequency components. 
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Clearly, the assumption that only the constant portion of the spectrum of the 
LFM pulse is significant is not valid. Inspection of Figs. 5.4a through 5.4d justifies our 
scheme to calculate the number of frequency components required to represent a LFM pulse 
shape given by Eq. (5.28). Note that Eq. (5.28) was derived for the case of rectangular 
amplitude modulated pulses. Since the rectangular-shaped function contains a severe 
discontinuity', we expect that the number of frequency components given by Eq. (5.28) will 
coincide to a worst case situation. Consequently, we would expect that if we chose a 
modulating function that did not have as large a discontinuity at the leading and trailing 
edges of the pulse, we would not require as many frequency components to accurately 
represent the signal. 

F. ALTERNATIVE AMPLITUDE MODULATION FUNCTIONS 

As discussed earlier in this chapter, although w-e are primarily interested in the 
rectangular amplitude modulating function, it is not the only one that will be considered. 
Of interest is the Hanning weighting function given by 



r A [0.5 + O.5cos(27tt / T p )] , |t|<T p /2 

a(t)=l 0 , !■[> Tp/ 2 



(5.32) 



and the Hamming w eighting function given by 

r A[0.54+0.46cos(27tt/T p )] , |t|<T p /2 
a(t)_1 0 , |t|> T p/ 2 



(5.33) 



where A is a constant and Tp is the pulse length (in seconds). It should be noted from 
Eqs. (5.32) and (5.33) that the pulse is assumed to be centered at t = 0, as was the case for 
the rectangular amplitude modulating function given by Eq. (5.16). Note that Eqs. (5.32) 
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and (5.33) are the actual pulse amplitude modulation functions and should not be confused 
with the concept of windowing functions applied to a block of data, i.e., the modulation 
functions are applied specifically to the pulse (of length Tp seconds) and not to the entire 
data record (of length T 0 seconds). 

Inspection of Eqs. (5.32) and (5.33) reveals that they are both a maximum at the center 
of the pulse, and decrease towards zero at the leading and trailing edge of the pulse in a 
"smooth" fashion. It can also be readily seen that these two equations are very nearly 
identical. Why then are we interested in two very similar amplitude modulation functions 
of this form? From the discussion in the previous section, we expect that it will require 
less frequency components to accurately represent the transmitted signal than was the case 
when a rectangular amplitude modulation function was used. To test this hypothesis, 
examples of amplitude modulated CW pulses are plotted as follows: 

* Fig. 5.5a. Hanning amplitude modulation, 50% duty cycle, 

• Fig. 5.5b. Hanning amplitude modulation, 8% duty cycle, 

• Fig. 5.6a. Hamming amplitude modulation, 50% duty cycle, 

* Fig. 5.6b. Hamming amplitude modulation, 8% duty cycle. 

In Figs. 5.5a and 5.6a, the pulse reconstructed from the calculated Fourier series 
coefficients of its complex envelope, denoted x ce (n), is compared with the actual sampled 
bandpass signal, denoted x(n), where x(n) is plotted as a dashed curve while x ce (n) is 
plotted as a solid curve. For plotting purposes, a third-degree polynomial curve fitting 
routine is used to interpolate between the data points. In Figs. 5.5b and 5.6b, just Xc e (n) 
is plotted to prevent the plot from becoming too cluttered. 
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Comparing Figs. 5.5a and 5.6a with Fig. 5.1a, it can be seen that indeed fewer 
frequency components are required to accurately represent the transmitted signal when the 
amplitude modulation function is of the form given by Eqs. (5.32) or (5.33), than if it is of 
the form given by Eq. (5. 16). Whereas we introduced the concept of Lanczos smoothing 
for the rectangular modulated pulse in an attempt to obtain a ''cleaner" representation of the 
signal, we find that this is totally unnecessary for the Hanning and Hamming modulation 
functions. The observations that can be made when comparing Figs. 5.5b and 5.6b with 
Fig. 5.3a are consistent with those just described. From these comparisons we can make 
the statement that the Hanning and Hamming modulated CW pulses require roughly half 
the number of frequency components for accurate representation than does the rectangular 
modulated pulse. No comparison of modulated LFM pulses was carried out as it was 
judged past the scope of interest for this thesis. 

Further comparison of Figs. 5.5 and 5.6 shows that for a given number of frequency 
components the representation of the Hamming amplitude modulated pulse is better than the 
Hanning modulated pulse. This means that there is closer agreement between the pulse 
reconstructed from the calculated Fourier series coefficients of its complex envelope and the 
actual sampled bandpass signal when the Hamming amplitude modulation function is used. 
In particular, the Hamming modulated pulse displays significantly less ringing before/after 
the leading/trailing edges of the pulse for the lower duty cycle case shown in Figs. 5.5b 
and 5.6b. 

The advantage in reducing the number of frequency components required to represent a 
signal is a simple one. Since the computer simulation must compute a value of the overall 
system complex frequency response for each frequency of interest, and this operation 
comprises the majority of the computational effort, a decrease in the number of required 
frequencies will mean an almost proportional drop in computation time, i.e., if the 
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fl: 1.0 IP- 12.5 MSEC PRF: 10.0 HZ FC: 100.0 HZ 

Fig. 5.5a CW pulse, 50% duty cycle, Hanning amplitude modulation. 




CF1SE: FIG 5.5B WFlVEr URM •■CW PULSE ( NFREO: 39 ) 

0: 1.0 TP: 10.0 MSEC FRF: 2.0 HZ FC: 100.0 HZ 



Fig. 5.5b CW pulse, 8% duty cycle, Hanning amplitude modulation. 
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TRANSMITTED ELECTRICnL SIGNAL X 




Fig. 5.6a C\V pulse, 50% duty cycle, Hamming amplitude modulation. 




R: 1.0 TF: -!0.0 MSEC FP^: 2.0 HZ FC: 40C.O HZ 

Fig. 5.6b C\Y pulse, 8% duty cycle, Hamming amplitude modulation. 
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number of frequencies is reduced by half, then the program run time will be approximately 
reduced by half. 

However, this is not the main reason that we consider amplitude modulation functions 
of this form. As stated earlier, we are primarily interested in the rectangular amplitude 
modulating function. Our interest in the Hamming and Hanning functions stems from the 
fact that many currently available computer models dealing with sound propagation in the 
ocean deal with transmitted waveforms having these particular amplitude modulation 
functions, presumably in an attempt to reduce computation times. If we then also have the 
option of using these same modulation functions, we will be in a better position to compare 
our results and performance with other ocean medium propagation models, e.g., SAFARI 
[Ref. 7] . 
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VI. COMPUTER SIMULATION RESULTS FOR THE 
ISOSPEED MEDIUM CASE 



A. THE OVERALL SYSTEM COMPLEX FREQUENCY RESPONSE FOR AN 
ISOSPEED MEDIUM 

The first step in validating the developments carried out in Chapters 2 through 4 is to 
apply the analysis to the simplest scenario, i.e., the free-space, isospeed medium case. 
Since the preceding development of the overall system complex frequency response given 
by Eq. (3.14) assumed a free-space problem, all we need to do is apply the condition of an 
isospeed medium to the ocean medium transfer function given by Eq. (4.34). As 
previously discussed, the phase integral of the ocean medium transfer function has a closed 
form expression in the isospeed case, which results in the relationship 



= -M k Y(yp)| r -j^Y(yoXy-y 0 )pj27tfY(yo)(y-yn) 
2 k Y(y 0 ) -/| k y(y)| 



where 



k Y (y 0 ) = ±2^[(f/c o ) 2 -f?] , f?£(f/c 0 ) 2 (6.2a) 

2 7 ^ 2 2 

k Y (y 0 ) = + j2rc[f r -(f/c 0 n , f r >(f/c Q ), (6.2b) 

k Y (y 0 ) = 27tf Y (y 0 ) = 27tf Y , and c 0 is the constant speed of sound in the isospeed 

medium. Equation (6.1), then, is our ocean medium transfer function for a free-space, 
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isospeed medium. Note that the ocean medium transfer function given by Eq. (6.1) 
contains no integral operations due to there being a closed form expression for the phase 
integral of Eq. (4.34). This result is significant since numerical integration routines are 
very time consuming. Therefore, we expect that the computational effort required for the 
isospeed case will form a baseline against which to compare other cases that do not have 
closed form solutions. 

It can be easily seen from Eqs. (6.2a) and (6.2b) that for an isospeed medium, the 
propagation vectors in the Y directions, k Y (yo ) and k y (y), are equal. The ray acoustics' 
interpretation of this is that rays travel in straight lines and will not bend. It should be 
remembered from the development in Chapter 3, that Eq. (6.2a) represents propagating 
waves while Eq. (6.2b) represents evanescent waves. 

Returning to our expression for the overall system complex frequency response given 
by 



we notice that the region of integration with respect to f r is from zero to infinity. Equation 
(3.14) is referred to as a full-wave solution for our ocean medium propagation problem, 
and can be thought of as an infinite summation of ray acoustic solutions [Ref. 7:p. 36]. 
From Eqs. (6.2a) and (6.2b), it can be seen that this region includes both propagating and 
evanescent waves. Since the propagation vector in the Y direction has different forms 
depending on whether the wave is propagating or not, we will split Eq.(4.33) into two 
integrals. One integral will cover the region of propagating waves, i.e., integrate with 



H(f,r) = 2 jc 




c mn(f) (f,f r ,yo;y) J()( 27 tf r r m ) 






(3.14) 
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respect to f r from zero to (f / c 0 ), while the second will cover the region of evanescent 
waves, i.e., integrate with respect to f r from (f / c 0 ) to infinity. 

As discussed in Chapter 4, our general solution for the case of a unit amplitude point 
source in an unbounded, homogeneous medium should collapse into the form of the free- 
space Green's function for an isospeed medium given by 



<p(r) = 



4tiR 



(4.25) 



In this chapter, we will refer to Eq. (4.25) as the exact solution for the free-space, isospeed 
medium since no simplifying approximations were made to obtain this result. Examination 
of Eq. (4.25) shows that under these circumstances, the propagating wave undergoes 
spherical spreading, i.e., the amplitude of the acoustic field is inversely proportional to 
radial distance between the transmitter and the field point. This radial distance is 
commonly referred to as the Range-Line-Of-Sight (RLOS). 

Reviewing the development leading up to the expression for the overall system complex 
frequency response due to a unit amplitude point source given by 

H(f,r) = 2jt f H M (f,f r ,y o ;y)J 0 (2nf r r m )e' j2 ” fY(y<>)<y ' y " ) f r df r , (4.33) 

•'o 

remember that it is mathematically the same as the acoustic field given by Eq. (4.25). As a 
result, the overall system complex frequency response should behave the same as the 
acoustic field due to a point source. Accordingly, if we plot RLOS versus the magnitude 
of the overall system complex frequency response on log-log axes for a particular 
frequency component of the transmitted signal, we would expect the result to be a straight 
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line. Further, due to the isospeed nature of the medium, we would expect all propagating 
frequency components to behave the same since we are only concerned in this simulation 
with the effects of spreading and have ignored losses due to absorption. 

B. RESULTS OBTAINED BY IGNORING EVANESCENT WAVES 

An initial investigation was carried out for an isospeed medium with a sound speed of 
1475 m/sec. For a transmit signal, we used the waveform of Fig. 5.1b, which was a 
Lanczos smoothed, rectangular amplitude weighted CW pulse with the following 
parameters : 

• AMP = 1.0, TP = 12.5 msec, PRF = 40.0 Hz, FC = 400.0 Hz 

• minimum frequency = 160.0 Hz, maximum frequency = 640 Hz. 

For this scenario, the longest wavelength (which corresponds to the lowest frequency 
component present) is equal to 9.219 m, while the shortest wavelength (which corresponds 
to the highest frequency component present) is equal to 2.305 m. Since the evanescent 
waves are exponentially decaying instead of being oscillatory, they are considered as being 
non-propagating and are usually ignored at distances greater than a couple of wavelengths 
from their source. For the case just described, we would expect evanescent waves to be 
insignificant at ranges greater than say fifty meters. Applying this assumption to our 
expression for the overall system complex frequency response, i.e., ignoring the 
integration over the region of evanescent waves, results in 

f iUo 

H(f.r)-2* I H M (f,f r .y o ;y)J 0 (27tf r r m )e' j2 ’ tfY °' oXy ' y ” ) f I df r , (6.3) 
*'0 
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which is an approximation of the full-wave solution since the region of integration is no 
longer infinite. We will refer to any overall system complex frequency response that has 
such a limited region of integration as our approximate (full-wave) solution to the ocean 
medium propagation problem. 

It should be stressed at this point that the computer implementation of the overall 
system complex frequency response has been kept as general as possible. This means that 
instead of having to program separate forms of Eq. (3.14) for different scenarios, e.g., Eq. 
(6.3), just the generalized form of the overall system complex frequency response given by 
Eq. (3.14) is implemented in the computer code. This generalized form can then be 
modified as required, by a number of user defined variables. To obtain the overall system 
complex frequency response given by Eq. (6.3), we merely set the user defined variables 

• transmit array : MT = 1, NT = 1, STEER = FALSE 

• integration routine : RATIO = 1.0 

• sound speed profile : NGRAD = 0 

where MT and NT describe the dimensions (number of elements) of the transmit array, as 
discussed in Chapter 2, and STEER is a logical variable that designates whether any 
beamsteering is to be done. If no beamsteering is to be done, then the transmit array 
complex weights, given by Eq. (2.25), are set equal to unity. The variable, RATIO, 
defined as 



RATIO = f r c(y 0 ) / f , 



(6.4) 
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determines the upper limit of the integration with respect to f r , while the variable NGRAD 
set equal to zero designates an isospeed medium and the program uses the closed form 
expression for the phase integral appearing in Eq. (3.14). 

To demonstrate the validity of the preceding discussion, we need to determine if our 
derived overall system complex frequency response does in fact reduce to the isospeed, 
free-space, Green's function and show the expected (1 / RLOS) fall-off in the magnitude of 
the acoustic field. In other words, we need to compare the results given by the 
approximate full- wave and exact solutions. Accordingly, the magnitude of the overall 
system complex frequency response versus range-line-of-sight is plotted in Fig. 6.1 for 
three different frequencies. The frequencies chosen correspond to the minimum, carrier, 
and maximum frequencies of the transmitted signal. Additionally, the approximate full- 
wave and exact values of the magnitude of the acoustic field at the carrier frequency of 400 
Hz and for the values of RLOS used in Fig. 6.1 are detailed later in Table 6.3. 

As can be seen from Fig. 6. 1 , the magnitude of the acoustic field has been determined 
for values of RLOS ranging from ten meters to ten thousand meters. This is equivalent to 
a range of one to thousands of wavelengths, depending on which of the frequency 
components we are considering. From our preceding discussion, we would expect that 
the acoustic field calculated using the approximate full-wave solution, given by Eq. (6.3), 
will not reduce exactly to the isospeed, free-space, Green's function for ranges less than 
about fifty meters. This is because we have assumed that evanescent waves have a 
noticeable effect in this region, which Eq. (6.3) ignores. Conversely, would we expect 
close agreement for ranges greater than fifty meters as we have assumed that the effect of 
evanescent waves in this region is negligible. In terms of the log-log plot of magnitude 
versus RLOS shown in Fig. 6.1, we would expect the graph to be a straight line for ranges 
greater than fifty meters, and to deviate away from this ideal for shorter ranges. 
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Acoustic Field 

Velocity Potential ( m**2 / sec ) 




Fig. 6.1 Magnitude of the acoustic field versus RLOS when 
the effects of evanescent waves are ignored. 
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It can readily be observed from Fig. 6.1 that for long ranges, the graph does indeed 
approximate a straight line with the expected slope, and for short ranges, it deviates from 
this straight line. These results appear, at first glance, to be what we expected. However, 
note the range at which the plotted curve starts to approximate a straight line. This range is 
about five hundred meters, or ten times the range that we thought would be significant. 
Alternatively, this range represents between 54 and 216 wavelengths, depending on which 
of the three frequency components we are considering. This would suggest that the effects 
of evanescent waves are indeed still significant at ranges greater than a couple of 
wavelengths from their source. 

Consider also, the actual values of the approximate acoustic field that are contained in 
Table 6.3. Even at a range of ten thousand meters, the approximate full-wave and exact 
solutions calculated for the carrier frequency, still differ by almost 1%. For the carrier 
frequency, this range corresponds to 2,712 wavelengths, which is well beyond where we 
expected evanescent waves to have an effect. 

Close inspection of Fig. 6.1 does indicate, however, that the low frequency term, i.e., 
160 Hz, is the one that is behaving the most erratically. This is what we would expect 
with regards to the deviation from the straight line predicted by the exact solution since this 
frequency corresponds to the longest wavelength. Consequently, for any given distance, 
this is the evanescent component that would have decayed the least. Since we are ignoring 
evanescent waves in this case, this is the frequency component that should deviate most 
from the exact result. Based on these observations and trends, we will assume that our 
derivation of the overall system complex frequency response has not been proved in error 
so far, and that it is the assumption regarding the significance of evanescent waves that is 
incorrect. 
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Even though our approximate solution results deviate from those predicted by the exact 
solution, is the amount of deviation really significant ? Remember that one of the outputs 
that we wanted from the computer simulation was the time-domain electrical signal at each 
element in the receive array, which would be used to show pulse distortion due to the 
effects of dispersion, and as input data to space-time signal processing algorithms. To 
demonstrate this pulse distortion, we need to compare the original transmitted electrical 
signal with the received electrical signal and observe the changes that propagation through 
the ocean medium has caused. Such a comparison can be made using Figs. 6.2a and 
6.2b. Note that 6.2b is the received electrical signal at a RLOS of 100 m. Referring to 
Fig. 6.1, we see that this value of RLOS is one for which the calculated acoustic field 
varies significantly with frequency. This particular RLOS was chosen to represent a near 
worst case situation and will demonstrate the significance, if any, of the inexactness of the 
approximate full-wave solution when evanescent waves are ignored. 

It can be seen from comparing Figs. 6.2a and 6.2b, that the output electrical signal is 
indeed significantly different from the transmitted signal. Recall that there should be no 
dispersion in an isospeed, free-space medium. This means that the received electrical 
signal should be a scaled replica of the transmitted electrical signal, i.e., identical in pulse 
shape. The problem with the approximate full-wave solution that produced the output 
electrical signal shown in Fig. 6.2b is that, if we repeated the problem for a non-isospeed 
medium in which we do expect pulse distortion due to dispersion, we would not know 
how much of the distortion was due to the dispersion and how much was due to the 
approximate nature of our solution. 

Next, consider the previous problem repeated for a RLOS of 1,000 m. From Fig. 6.1 
it can be seen that for this range, the acoustic field behaves a lot more as predicted and is 
less dependent on frequency. This means that there is much closer agreement between the 
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R: 1.0 TF: 12.5 MSEC FRF: 10.0 HZ FC: 100.0 HZ 

Fig. 6.2a Transmitted electrical signal. 




YT: 1000.0 M XR: 0.0 11 YR: 1050.0 H ZR: 86.5 M 
RLOS: 100.0 M BLOS: 60.0 DEG 



Fig. 6.2b Received electrical signal, RLOS = 100 m, when the 
effects of evanescent waves are ignored. 
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approximate full-wave and exact solutions than existed for RLOS = 100 m. The output 
electrical signal for this case is shown in Fig 6.2c, where the transmitted electrical signal is 
the same as for the previous case shown in Fig. 6.2b. As expected, the output electrical 
signal indeed looks a lot more like the transmitted electrical signal than was the case in Fig. 
6.2a. In fact the amount of distortion can be considered imperceptible for RLOS = 1,000 
m, and, for all practical purposes, Fig. 6.2c is a scaled version of the transmitted signal. 

Based on these results for the received time-domain electrical signal, it appears that the 
approximate full-wave solution is adequate for demonstrating the effects of pulse distortion 
for signals containing no frequency components below 160 Ilz, where the receive element 
is at a range greater than 1,(XX) m from the source. This situation should be acceptable for 
long range propagation problems, but for any short range calculations, it is obviously 
inadequate. 




CRSEsHOM CW 1.0 

YT: 1000.0 M XR : 0.0 H YR: 1500.0 M ZR: 666.0 M 
RL05: 1000.0 M BLOS: 60.0 050 



Fig. 6.2c Received electrical signal, RLOS = 1,000 m, when 
the effects of evanescent waves are ignored. 
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Consider next, the other use of the time-domain electrical signal at each element in the 
receive array, that is, as input data to space-time signal processing algorithms. As stated 
earlier, the computer code developed for this ocean medium propagation problem does not 
perform any signal processing on the received electrical signal at each point in the array, but 
rather, the data is written to a file so that it can be operated on by ancillary programs. As a 
further validation of the generated results, we can operate on the time-domain data with a 
beamforming algorithm that localizes the position of the source, and compare the computed 
and actual source position. The algorithm used was an adaptive beamforming and non- 
linear estimation algorithm developed by Breckenridge [Ref. 16] that localizes a source in 
range, depression angle, and azimuthal angle, according to the geometry outlined in Fig. 
6.3. From Fig. 6.3 it should be noted that the localization of the source, calculated by this 
adaptive algorithm, is relative to the center element in the receive array. 

The data generated from the two cases, RLOS = 100 m and RLOS = 1,000 m, was 
processed with this algorithm using the output electrical signal from each of the elements in 
a 3x3 receive array, and the results are shown in Table 6.1. The values RLOS (X) and 
RLOS (Y) detailed in the table correspond to estimates of the RLOS calculated from 
information along the X and Y axes, respectively. From the results contained in Table 
6.1, we can see that the estimated values from the beamforming algorithm are closer to the 
true values for the case RLOS = 1,000 m than for the case RLOS = 100 m, which is 
consistent with our previous observations in this section. 

An interesting feature of the results is the way that RLOS (Y) differs radically from the 
true value of RLOS, while RLOS (X) is fairly close to the true value. Note that this occurs 
for both values of RLOS considered, that is, RLOS (Y) does not improve as the value of 
RLOS is increased. This phenomena is in contrast to the trends observed earlier in this 
section. The fact that the effects of evanescent waves are being ignored means that we are. 
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Fig. 6.3 Geometry for adaptive beamforming algorithm. 
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Table 6.1 Adaptive beamforming results when 
evanescent waves are ignored. 
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in essence affecting the wave solution in the Y direction to a much greater degree than in the 
X direction. Recall that our form of the overall system complex frequency response, given 
by Eq. (3.14), was developed assuming a sound-speed profile that was a function of 
depth. In our development we showed that this meant Eq. (3.14) was axisymmetric. If 
the overall system complex frequency response is axisymmetric, then the relative behavior 
of the acoustic field in the X direction should be minimally affected by the behavior of the 
acoustic field in the Y direction. Consequently, if there is a problem with the description 
of the full-wave solution in the Y direction, as in this case, the value of RLOS (X) should 
still be fairly accurate, which is the result that we observed. 

These results also demonstrate that while the time-domain output electrical signal from a 
single receive array element may look correct to the eye and suggest that the accuracy of the 
generated data is adequate, it may in fact not be accurate or correct enough for signal 
processing purposes. This suggests that observing the output electrical signal is not a 
good way of validating the performance of the approximate full-wave solution. 

The remaining output that we wanted from this computer simulation was the magnitude 
and phase of the complex acoustic field incident upon a receive array as a function of 
frequency and spatial coordinates. A convenient way of viewing these results, apart from 
a three-dimensional plot, is to tabulate the magnitude and phase of the acoustic field in the 
same order as the elements of the receive array. Consider a 3x3 element receive array. 
We will adopt the element designation scheme shown in Fig. 6.4 to number the elements 
( m,n ) contained in the array. 

The values of the magnitude and phase (in degrees) of the acoustic field are shown 
in Tables 6.2a and 6.2b for the cases RLOS = 100 m and RLOS = 1,000 m, respectively. 
Note that the first value is the magnitude and the second value is the phase. Given the 
placement of the receive array relative to the source, i.e., XR = 0 and YR = 1025 m, which 
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Fig. 6.4 Element designation for a 3x3 receive array. 

is below the source, we would expect the magnitude of the acoustic Field to decrease with 
increasing n. This corresponds to reading down the table, or equivalently, increasing the 
depth of the receive element or field point. We would also expect the value of the acoustic 
field to decrease with increasing absolute value of m. This corresponds to reading across 
the table out from the center, or equivalently, increasing the displacement in the X direction 
of the element or field point. We expect the magnitude of the acoustic field to decrease in 
both of these situations since their effect is to increase the radial range between the source 
and the receive element , or field point. 

From these two tables it can be seen that all the values for the acoustic field do not 
behave as expected. In Table 6.2a, the first two rows of results have the magnitude of the 
acoustic field increasing as we move away from m = 0. In Table 6.2b on the other hand, it 
is the third or final row of values that behaves in this fashion. In both cases however, the 
magnitude of the acoustic field decreases as the depth is increased, as expected. As for the 



93 
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-1 


8.69822 e-4 
-168.608 


8.68689 e-4 
-167.994 


8.69822 e-4 
-168.608 
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8.47479 e-4 
139.860 


8.47246 e-4 
140.423 


8.47479 e-4 
139.860 


1 


7.82107 e-4 
84.800 


7.82997 e-4 
85.376 


7.82107 e-4 
84.800 


Table 6.2a Acoustic field at array elements, RLOS = 100 


n/m 


, 


0 




-1 


8.00331 e-5 
169.065 


8.00375 e-5 
169.130 


8.00331 e-5 
169.065 


0 


7.97398 e-5 
113.082 


7.97420 e-5 
113.150 


7.97398 e-5 
113.082 


1 


7.93029 e-5 
56.758 


7.93010 e-5 
56.826 


7.93029 e-5 
56.758 



Table 6.2b Acoustic field at array elements, RLOS = 1,000 m. 
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results obtained from the use of the adaptive beamforming algorithm, observing the values 
of the complex acoustic field as a function of frequency and spatial coordinates alerts us to 
the fact that the data produced by our approximate full-wave solution is inaccurate. 

From the results of this section, we can formulate two primary conclusions. First, the 
assumption that the effects of evanescent waves past a distance of a couple wavelengths is 
insignificant is an invalid one, at least for the situation we are modelling and for the full- 
wave method used. Second, comparing the transmitted and received electrical signals for 
the free-space, isospeed medium case is not sufficient to determine whether or not the 
approximate full-wave solution, is accurate. To gauge the performance of the approximate 
solution we need to either observe the acoustic field as a function of spatial coordinates for 
an appreciation for whether the solution is behaving correctly, or apply signal processing 
techniques to obtain some quantitative measure for the accuracy of the solution. 

C. RESULTS OBTAINED BY INCLUDING EVANESCENT WAVES 

To include the effects of evanescent waves, we need to extend the region of integration 
with respect to f r in Eq. (6.3). As was previously discussed, the region of evanescent 
waves extends from f r equal to (f / c G ) to infinity. But is it necessary to integrate over this 
entire region ? After some trial and error, it was found that integrating with respect to f r up 
to a limit of 1.2 (f / c 0 ) was adequate to give very good results down to a RLOS of 10 m. 
Recall that evanescent waves have more effect at short ranges. Therefore, if this upper limit 
of integration is sufficient for a RLOS of 10 m, then it will also be sufficient for longer 
ranges. Since no practical reason could be thought of as to why we would be interested in 
ranges less than 10 m, this is the upper limit of integration that we adopted and is that used 
to obtain the results of this section. Therefore, our approximate full-wave solution is now 
given by 
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f r df r . (6.5) 



Recall from our previous discussion that the region of integration in Eq. (6.5) is, for 
computational purposes, divided into two regions to account for the change in the form of 
the propagation vector in the Y direction, k y (y), when f r changes from being less than to 
greater than f / c 0 . This empirically determined value for the upper limit of integration 
agrees with observations made by Schmidt [Ref. 7:p. 31] with regard to wavenumber 
integrations. 

For the analysis contained in this section, the scenario is the same as for the previous 
case where we ignored the effects of evanescent waves, i.e., same transmit signal, array 
geometry and ocean medium. To allow for this extended region of integration, all that has 
to be done is to set the variable RATIO equal to 1.2 instead of 1.0. This demonstrates one 
of the advantages in having kept our initial analysis and computer code implementation as 
general as possible. It is a very simple matter to alter the problem that requires no 
programming changes. 

The magnitude of the overall system complex frequency response was determined for a 
range of values for RLOS and for three frequencies corresponding to the minimum, carrier, 
and maximum frequency components of the transmitted signal. It was found, however, 
that, unlike the previous case, there was no discemable frequency dependence in the values 
for the magnitude of the overall complex frequency response determined from the 
approximate full- wave solution given by Eq. (6.5). As was previously discussed, this 
frequency independence is the result predicted from the exact form of the isospeed, free- 
space Green's function. Since the results for the three frequency components were 
identical, only the results for the carrier frequency, 400 Hz, are shown in the log-log plot 
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of acoustic field versus RLOS shown in Fig. 6.5. Additionally, the approximate full-wave 
and exact values of the magnitude of the acoustic field for this frequency and values of 
RLOS used in Fig. 6.5 are detailed in Table 6.3. 

From Fig. 6.5, it can be observed that the log-log plot of the acoustic field magnitude 
versus RLOS is identical in slope to the straight line predicted by the isospeed, free-space. 
Green's function, over the entire range of values of RLOS considered. Comparing Figs. 
6.1 and 6.5, it can easily be observed that there has been a radical improvement due to the 
inclusion of a small region of evanescent waves. The results from Fig. 6.5 suggest that 
the upper limit of integration chosen in the approximate full-wave solution given by Eq. 
(6.5) is adequate for values of RLOS down to at least 10 m. 

Inspection of Table 6.3 shows that with the inclusion of this region of evanescent 
waves, the magnitude of the acoustic field given by the approximate full-wave solution in 
fact agrees with the exact values predicted by the isospeed, free-space, Green's function to 
within five significant figures. The comparison was limited to this number of significant 
figures since the numerical integration was only carried out to sufficient terms to give a 
maximum of a 0.01% relative error in the value of the integral. Alternatively stated, the 
value of the acoustic field calculated from Eq. (6.5) is accurate to within ±0.01%. 

Let us now consider the displaying of the received electrical signal for the purposes of 
showing pulse distortion due to dispersion. Comparing Figs. 6.6a and 6.6b, it can be 
observed that there is no discemable pulse distortion, which is what theory predicts for 
wave propagation in a free-space, isospeed medium. This result was also expected from 
experience with the analysis of the previous case, in which we ignored the effects of 
evanescent waves. In that case, even though there was some frequency dependence in the 
value of the acoustic field calculated using the approximate full- wave solution of Eq. (6.3), 
there was very little discemable distortion in the received electrical signal. For the acoustic 
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Acoustic Field 

Velocity Potential ( m**2 / sec ) 




Fig. 6.5 Magnitude of the Acoustic field versus RLOS when 
the effects of evanescent waves are included. 



98 



RLOS(m) 


l/4nR 


1 H(f,r) 1 


1 H(f,r) 1 






( no evanescent 
waves) 


( incl. evanescent 
waves ) 


10 


7.9577 e-3 


7.9702 e-3 


7.9577 e-3 


50 


1.5915 e-3 


1.3983 e-3 


1.5915 e-3 


100 


7.9577 e-4 


8.4725 e-4 


7.9577 e-4 


200 


3.9789 e-4 


3.9676 e-4 


3.9789 e-4 


300 


2.6526 e-4 


2.5923 e-4 


2.6526 e-4 


500 


1.5915 e-4 


1.5663 e-4 


1.5915 e-4 


1,000 


7.9577 e-5 


7.9742 e-5 


7.9577 e-5 


2,000 


3.9789 e-5 


3.9091 e-5 


3.9789 e-5 


3,000 


2.6526 e-5 


2.5978 e-5 


2.6526 e-5 


5,000 


1.5915 e-5 


1.5753 e-5 


1.5915 e-5 


10,000 


7.9577 e-6 


8.0182 e-6 


7.9577 e-6 



Table 6.3 Comparison between exact and predicted 
magnitude values of the acoustic field. 
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field calculated using the more accurate approximate full-wave solution given by Eq. (6.5), 
we would expect our received electrical signal to look even more like a scaled version of the 
transmit signal. This was, in fact, the observed result. 

The beamforming algorithm was then applied to the output electrical signals from a 3x3 
receive array, as in the previous case, with the results being summarized in Table 6.4. The 
first thing to notice about these results is that unlike the previous case that ignored 
evanescent waves, both RLOS (X) and RLOS (Y) give accurate values for the radial 
distance between the source and the center of the receive array. Additionally, the estimated 
position angles are very accurate. 

Finally, an inspection of the complex acoustic field at the elements of the receive array 
is shown in Table 6.5. It can be clearly seen that the values of the magnitude conform to 
the expected pattern, as described for the previous section. Note that the change in 
magnitude, reading down the table, is much greater than when reading across the table. 
While this may appear to be an error at first, it is in fact correct, and is due to the 
positioning of the array. Recall that the receive array is positioned such that XR = 0, 
YR = 1050 m, i.e., the receive array is 50 m below the transmit array, and RLOS = 100 m. 
Also recall that the distance between elements in the X and Y directions is the same. 
Accordingly, it can be easily shown that the change in radial distance between the source 
and a receive element is much larger for an incremental change in the Y direction than for 
the same change in the X direction. As a result, the change in the magnitude of the 
acoustic field will be similarly affected. 

From the results of this chapter so far, we have seen that it is not a valid assumption to 
ignore the effects of evanescent waves when developing an approximation of the full-wave 
solution to an ocean medium propagation problem. However, it appears that we only have 
to include a small portion of the evanescent region for an approximate full-wave solution to 
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RECEIVED ELECTRICAL SIGNAL Y 




A : 1.0 IP: 12.5 NSEC PRF: -10.0 HZ EC: 400.0 HZ 

Fig. 6.6a Transmitted electrical signal. 




CASE: HON CM 1.2 

YT: 1000.0 M VR: 0.0 N YR: 1050.0 N ZR: 65.6 M 
RL05: 100.0 M BLOS: 00.0 DEG 



Fig. 6.6b Received electrical signal, RLOS = 100 m, when the 
effects of evanescent waves are included! 
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be valid. From these results we will further assume that the region of integration with 
respect to f r , given by 



0 < f r < 1.2 (f/c 0 ), (6.6) 

for the overall system complex frequency response, given by Eq. (3.14), is adequate for 
sound-speed profiles other than the isospeed case investigated in this chapter. We will 
utilize this assumption in the following chapter, when we investigate the behavior of the 
acoustic field in a free-space medium characterized by a linear sound-speed profile with a 
single gradient. We have also seen that the accuracy of the approximate full-wave solution 
is a function of the numerical integration used to calculate the overall system complex 
frequency response. There is a trade-off between computation time and accuracy that will 
be addressed in the final section of this chapter. 

D. PLOTTING OF THE 3-D COMPLEX ACOUSTIC FIELD 

So far in this chapter, all our validation has taken place using a similar geometry, i.e., 
XR = 0 and YR > YT. This means that the receiver has always been deeper than the 
transmitter, and has had no offset in the X direction. To have any confidence that the 
theory and computer implementation of the overall system complex frequency response is 
correct, we need results for every permutation of the X and Y direction offsets of the 
receiver with respect to the source. In other words, we need to place the receiver in a 
variety of positions, e.g., positive X offset and above the source, to fully test the 
algorithm. 

Since one of the required outputs is the magnitude and phase of the complex acoustic 
field at the locations of the array elements for use by matched-field processing algorithms, a 
3-D plotting routine was implemented and added as a visual aid to help conceptualize what 



102 



RLOS 


RLOS (X) 


RLOS 00 


0 


0 


¥ 


W 


(m) 


(m) 


(m) 


(deg) 


(deg) 


(deg) 


(deg) 


true 


est 


est 


true 


est 


true 


est 


100 


99.984 


99.977 


30.000 


30.004 


270.00 


270.00 



Table 6.4 Adaptive beamforming results when 
evanescent waves are included. 
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8.00292 e-4 
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8.00346 e-4 
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8.00191 e-4 
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7.95722 e-4 
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7.95774 e-4 
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7.95722 e-4 
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7.91125 e-4 
79.909 


7.91177 e-4 
80.554 


7.91125 e-4 
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Table 6.5 Acoustic field at array elements, RLOS = 100 m. 
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the acoustic field looked like. This routine plots, on separate graphs, the magnitude and 
phase of the overall system complex frequency response for a particular frequency versus 
the relative position of each of the array elements. Recall from our earlier development that 
the value of the frequency spectrum of the acoustic field at any point is just the value of the 
overall system complex frequency response times the frequency spectrum of the transmitted 
electrical signal. By then looking at the plots, particularly those of the magnitude of the 
acoustic field, it can easily be determined if the field exhibits the appropriate orientation and 
curvature for the given relative positions of the source and the receive array. 

Note that the planar array of point sources we are calling the transmitter does not 
necessarily have to be a physical transducer array. It could represent sampled values of the 
acoustic field at some spatial location (or grid of locations) due to some other source. In 
this situation, the array could be thought of as detecting an acoustic signal and reradiating 
that acoustic signal without alteration. 

As was previously discussed, we need to generate a series of plots of the complex 
acoustic field for various combinations of source / receiver locations. A series of plots 
were generated for the six different permutations of the X and Y offsets of the receive array 
relative to the source. For ail these plots, the source was located at a depth of 1,000 m. 
Recall that the transmit array (or single point source in this case) is centered within the 
coordinate axis system such that XT and ZT are both equal to zero and YT is the depth 
below the sea surface (see Fig. 2.3). The receiver was a 9 x 9 planar array of point 
sources, which resulted in the computation of the overall system complex frequency 
response for 81 different spatial locations. These values of the overall system complex 
frequency response were calculated for a frequency of 400 Hz, which is the carrier 
frequency of the transmitted electrical signal shown in Figs. 6.2a and 6.6a. The RLOS, 
i.e., the radial distance between the source point and the center of the receive array, was 
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assumed to be only 50 m in order to reduce computer processing time. This short range 
was chosen so that we could process data a reasonable sized array which would allow us 
to observe curvature in the acoustic field, i.e., so that the plots would not just look like flat 
planes. The array size was chosen as being adequate to produce a 3-D plot, for the RLOS 
as previously chosen, in which the orientation and curvature of the acoustic field values 
was easily discemable. A constant speed of sound of 1475 m / sec was assumed for all 
cases. In the following 3-D plots, the X and Y offsets AY = YR - YT and AX = XR are 
used to describe the relative positions of the centers of the transmit and receive arrays. 

A series of six sets of plots were generated, showing both the magnitude and phase of 
the overall system complex frequency response, for the following positions of the receive 
array : 



Fig. 6.7 YR = 1000 m, XR = 0 m, (AY = 0, AX = 0) 

Fig. 6.8 YR = 1050 m, XR = 0 m, (AY = 50, AX = 0) 

Fig. 6.9 YR = 1010 m, XR = 20 m, (AY = 10, AX = 20) 

Fig. 6.10 YR = 980 m, XR = 10 m, (AY = - 20, AX = 10) 
Fig. 6.1 1 YR = 1020 m, XR = - 20 m, (AY = 20, AX = - 20) 
Fig. 6.12 YR = 990 m, XR = - 20 m, (AY = - 10, AX = - 20) 



Recall that the parameters, (XR, YR, ZR), describe the position of the center element of the 
receive array. 

Inspection of the magnitude plots of these figures clearly show that the computer code 
does indeed handle the different permutations of the X and Y offsets. Consider the 
magnitude plot shown in Fig. 6.10a. Since AY is negative and AX is positive, then from 
the perspective of the receiver looking towards the source, the receiver is located above and 
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CR5E = HOflOG CW FREQUENCY: 400.0 HZ 

YT: 1000.0 fl XR: 0.0 n YR: lOOO.O II ZR: 50.0 M 

DXR^ 1.152 M DYR: 1.152 fl 



Fig. 6.7a Magnitude of the overall system complex frequency response. 
(AY = 0 , AX = 0) 
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CRSE.-HOMOG CW FREQUENCY : 100. □ HZ 

YT: 1000.0 M XR: 0.0 M YR: 1000.0 M ZR: 50.0 M 

DXRi ■ 1.152 M DYR: 1.152 II 



Fig. 6.7b Phase of the overall system complex frequency response. 
(AY = 0, AX = 0) 
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CFISE : HOtlOG CW FREQUENCY: 400.0 HZ 

YT: 1000.0 M XR: 0.0 M YR: 1050.0 11 ZR: 86.6 M 

DXR: ' 1 . 152 M DYR: 1.152 11 



Fig. 6.8a Magnitude of the overall system complex frequency response. 
(AY = 50, AX = 0) 
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CRSE: HOriOG CW FREQUENCY: 400.0 HZ 

YT: 1000.0 M XR: 0.0 M YR: 1050.0 M ZR: 86.6 M 

DXR: 1.152 M DYR: 1.152 M 



Fig. 6.8b Phase of the overall system complex frequency response. 
(AY = 50, AX = 0) 
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CASE : HOMOG GW FREQUENCY: 400.0 HZ 

YT: 1000.0 H XR: 20.0 M YR: 1010.0 M ZR: 44.7 11 

DXR: '1.152 M DYR: 1.152 11 



Fig. 6.9a Magnitude of the overall system complex frequency response. 
(AY = 10, AX = 20) 
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CR5E: HOMOG CM 
YT: 1000.0 M 
DXR:’ 1.152 M 



FREQUENCY: 400.0 HZ 
XR: 20.0 H YR: 1010.0 M 
DYR: 1.152 fl 



ZR: 44.7 M 



Fig. 6.9b Phase of the overall system complex frequency response. 
(AY = 10, AX = 20) 
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CnSE:HOMOG CN 
YT: 1000.0 M 
DXRi 1.152 M 



FREQUENCY: 400.0 HZ 
XR: 10.0 M YR: 980.0 M 
DYR: 1.152 M 



ZR: 44.7 M 



Fig. 6.10a Magnitude of the overall system complex frequency response. 
(AY = - 20, AX = 10) 
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1 . 



CASE: HOMOG CW FREQUENCY: 4Q0.0 HZ 

YT: 1000.0 f1 XR: 10.0 M YR: 980.0 h ZR: 44.7 U 

DXR: 1.152 M DYR: 1.152 11 



Fig. 6.10b Phase of the overall system complex frequency response. 
(AY = -20, AX = 10) 
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CflSE:HOMOG CW FREQUENCY : 400.0 HZ 

YT: 1000.0 M XR: -20.0 M YR: 1020.0 M ZR: 41.2 M 

DXR: ' 1.152 M OYR: 1 . 152 H 



Fig. 6.1 la Magnitude of the overall system complex frequency response. 
(AY = 20, AX = - 20) 
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CR5E = HOhOG CM FREQUENCY: 100.0 HZ 
YT: 1000.0 II XR:-20.0 M YR: 1020.0 11 
DXR: ' 1.152 II DYR: 1.152 11 



ZR: 11.2 ri 



Fig. 6.1 lb Phase of the overall system complex frequency response. 
(AY = 20, AX = - 20) 
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CR5E: HOMOG CW FREQUENCY: 400.0 HZ 

YT: 1000.0 M XR:-20.0 M YR: 990.0 M ZR: 44.7 M 

DXR:' 1.152 M DYR: 1.152 M 



Fig. 6.12a Magnitude of the overall system complex frequency response. 
(AY = - 10, AX = - 20) 
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1 . 



CR5E : HOMOG CW FREQUENCY: 400.0 HZ 
YT: 1000.0 M XR:-20.0 II YR: 990.0 M 
DXRi ' 1.152 M DYR: 1.152 N 



ZR: 44.7 M 



Fig. 6.12b Phase of the overall system complex frequency response. 
(AY = -10, AX =-20) 
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to the left of the source. Remember also that this geometry relates to the center of the 
receive array. Given this geometry, we expect that if we move left of center, which 
corresponds to moving in a positive direction along the X axis or increasing the element 
index m (see Fig. 6.4), the smaller the acoustic Field would become. Conversely, if we 
move to the right of center, the radial distance between the source and the field point 
decreases and the acoustic field increases. Also given this geometry, as we go deeper, 
which corresponds to movement in a positive direction along the Y axis or increasing the 
array element index n, the radial distance between the source and the field point decreases 
and the acoustic field increases. This behavior is clearly demonstrated in Fig. 6.10a. 
Note that in this explanation, it was assumed that the field points did not cross the dividing 
lines where either AX or AY is zero. 

Note the magnitude plot of Fig. 6.7a, and the way the surface caves in by 
approximately 10% around the center element along the Y axis, i.e., n = 0. This 
corresponds to the center of the receiver being at the same depth as the source, i.e., axis- 
axis propagation. This phenomenon is caused by the fact that our expression for the ocean 
medium transfer function is based on the WKB approximation and is not exact. It was 
observed in Chapter 4 that as the propagation vector in the Y direction ky (y) approached 
zero, the approximate WKB solution is not valid since the solution approaches infinity. 
This can also be thought of in terms of the simplifying assumption given by Eq. (4.17) 
being invalid in this region. 

Inspection of Eq. (6.5) reveals that the integration with respect to f r will cause the 
ocean medium transfer function to be calculated for a wide range of values of the 
propagation vector in the Y direction for all spatial locations. Why is it only significant 
when the transmitter and receiver are at close to the same depth? By thinking of the 
problem from the perspective of ray acoustics, with the full- wave solution consisting of an 
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infinite summation of rays, it can be seen that as the transmitter and receiver get closer 
together in depth, the eigenray which connects the two gets closer to being horizontal and 
hence, closer to the invalid region of the WKB approximation. Since the eigenray and, 
those close are a major component of the acoustic signal at the receiver, it is when they get 
close to the invalid region that the effect becomes noticeable. 

This phenomenon also demonstrates an important point about ray acoustics. In ray 
acoustics, rays other than the eigenrays or those in very close proximity are often ignored 
and thought of as being unimportant. As a result, ray acoustics based on the WKB 
approximation would be of no use in an along axis, isospeed propagation problem since the 
resulting acoustic signal would be infinite. However, we have seen that the full-wave 
solution for an invalid region, by considering an infinite set of rays, was still able to give a 
solution to within 10% of the actual value. The full-wave solution is so called because it 
does not make any assumptions about what is important, and integrates over the entire 
region of propagating and evanescent acoustic wave contributions. We have shown that 
this infinite integration allows the solution to be equated to the exact solution for the free- 
space, isospeed problem and gives accurate values for the magnitude and phase of the 
complex acoustic field at a given field point. 

G. VALIDATION OF THE TRANSMIT ARRAY FAR-FIELD DIRECTIVITY 
FUNCTION 

Recall from Chapters 2 and 3 the discussion regarding the far-field directivity function, 
or beam pattern, of a planar array of point sources located in the XY plane. If our 
assumption about the axisymmetric nature of the far-field directivity function is correct, 
then applying the analysis of the preceding sections in this chapter to the overall system 
complex frequency response given by Eq. (3.14) yields 
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1.2 ( f / c„) mj. 



H(f,r) = 2 tc 



S c mn (f)H M (f,f r ,y 0 ;y)Jo(27cf r r m ) 



0 



m = -MT n = -NT 




(6.7) 



which should be a valid approximate full- wave solution for our ocean medium propagation 
problem. 

To validate Eq. (6.7), we will compute the value of the acoustic field at a number of 
field points (with a RLOS = 50 m) due to a 3x3 element transmit array and compare the 
results from exact calculations and the approximate full-wave solution. The transmitted 
electrical signal will be the Lanczos smoothed, rectangular weighted, CW pulse used 
extensively throughout this chapter. Recall that this signal has a maximum frequency 
component of 640 Hz and a carrier frequency of 400 Hz. For the purposes of this 
investigation, we will only consider the carrier frequency of the transmitted signal. 

First we need to check if we are in the far field. The intereiement spacing, in both the 
X and Y directions, of the transmit array is equal to one half of the wavelength 
corresponding to the highest frequency component of the transmit signal. Ziomek [Ref. 
3:pp. 118-120] shows that this value for the intereiement spacing is a necessary condition 
for avoiding grating lobes under all conditions of beam steering for planar arrays consisting 
of MT x NT (odd) elements. This is the default value used in the computer 
implementation if no user-defined intereiement spacing is specified, and is the value for 
which we will conduct this analysis. Combined with Eq. (2.23), this yields the condition 



RLOS > 




( 6 . 8 ) 
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which, if satisfied, places the field point in the far field of the transmit array and the far- 
field directivity function is the appropriate description of the beam pattern for this scenario. 
For this case, the right-hand side of Eq. (6.8) is equal to 2.263 m, which places the RLOS 
of 50 m well into the far field. 

If we further assume that the beam pattern is not steered and there is unit amplitude 
rectangular weighting of the array elements, it can be shown that Eqs. (2.28) and (2.30) 
combine to yield 



D . _ sin[u7t 15/16] sin[v7t 15/16] 
tU ' ~ sin[u7t5/16] sin[v;t 5/16] ’ 



( 6 . 8 ) 



where the direction cosines u and v are calculated from 



AX 
U RLOS 



(6.9a) 



and 



AY 

V RLOS * 



(6.10b) 



Values of the magnitude of the overall system complex frequency response for various 
combinations of AX and AY are detailed in Table 6.6. The exact value of the acoustic field 
is calculated by multiplying the relevant value of the far-field directivity function, given by 
Eq. (6.8), with the exact value of the acoustic field due to a single point source at the same 
RLOS calculated from Eq. (4.25). The simulated value of the overall system complex 
frequency response, on the other hand, is calculated using the approximate-full wave 
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solution given by Eq. (6.7). As can be easily observed from these results, the value of the 
acoustic field calculated using the approximate full-wave solution behaves well and is 
within 0.07% of the exact value for a wide range of direction cosine values. 

Recall from the discussion of the results detailed in Table 6.3 that the numerical 
integration of the overall system complex frequency response was carried out to sufficient 
terms to give a maximum of a 0.01% relative error. Note that the errors detailed in Table 
6.6 are larger than this value. Does this mean that there is an error with the assumption 
that the beam pattern is axisymmetric ? The answer is no. The increased error is a result 
of the manner in which Eq. (6.7) is numerically computed. To aid in the convergence of 
the numerical integration routine, the summations of Eq. (6.7) are taken outside of the 
integral. Therefore, each integration within the summation, with its associated error, 
contributes to the total error of the calculated overall system complex frequency response. 

Next, consider a steered beam pattern. The default value for the direction of 
beamsteering in the computer implementation is the line-of-sight direction between the 
center of the transmit array and the field point. In other words, the beam is steered 
towards the field point. This default value is used if no user-defined direction is specified, 
and this is the case we will investigate for this analysis. Since the beam is steered towards 
the field point, it can be shown that the far-field directivity function given by Eq. (2.28) in 
the direction of the field point reduces to a constant equal to the number of point sources in 
the anray. In other words, the array sources coherently add in the direction of the center of 
the main lobe of the beam pattern. As a result, the far-field directivity function of our 
steered transmit array has the value 9 in the direction of the field point. The values of the 
exact and approximate full-wave solutions of the acoustic field due to this steered array are 
detailed in Table 6.7. As can be observed, like the results contained in Table 6.6, the 
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exact and approximate results show close agreement, with errors being of the same order of 
magnitude. 

These results confirm that our analysis relating to the far-field directivity function and 
axisymmetric behavior of a planar array of point sources contained in the previous chapters 
was correct. The importance of this result is that it demonstrates that the assumption made 
in some other ocean propagation simulation programs, that axial symmetry can only be 
used for single point sources or vertical line arrays [Ref. 7] is incorrect 

F. IMPLEMENTATION OF NUMERICAL INTEGRATION 

To implement the integration of Eq. (3.14) with respect to f r , the IMSL10 integration 
routine, DQDAG, was used. This is a double precision numerical integration routine 
which uses a globally adaptive scheme based on Gauss-Konrad rules. Note that the 
integrand of Eq.(3.14) contains both a zero-order Bessel function of the first kind, 
implemented using the double precision IMSL10 function, DBSJO, and complex 
exponential terms. This means that the integrand is oscillatory, becoming more oscillatory 
as frequency and the range between the source and the field point are increased. As a 
consequence, a Gauss-Konrad rule is used with 30-61 points. 

The oscillatory nature of the integrand means that it is time consuming to compute 
numerically with a accuracy better than 0. 1 % - 0.5%. For our computer implementation 
we desired a relative error less than 0.01%, which was not achievable by numerically 
integrating between the appropriate limits, regardless of how much time was allowed. To 
achieve this error criterion, it was necessary to subdivide the original regions of integration 
into 20 subintervals, and to compute each interval separately and sum the result. 

It has been mentioned on numerous occasions that numerical integration is a time 
consuming operation. To illustrate this fact, some typical run times are detailed in Table 
6.8. Note that this program has been implemented in FORTRAN 77 and was run on the 
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AX 


Dr(f,fx) 


AY 


DrCf.fy) 


1 H(f,r) 1 
( exact ) 


1 H(f,r) 1 
( simulated ) 


% error 


0 


3.0000 


25 


2.1111 


1.0079 e-2 


1.0081 e-2 


0.020 


20 


2.4142 


-20 


2.4142 


9.2758 e-3 


9.2795 e-3 


0.040 


- 10 


2.8478 


10 


2.8478 


1.2907 e-2 


1.2903 e-2 


- 0.031 


-10 


2.8478 


-30 


1.7654 


8.0015 e-3 


8.0038 e-3 


0.029 


35 


1.3902 


10 


2.8478 


6.3008 e-3 


6.3050 e-3 


0.067 



Table 6.6 Acoustic Held resulting from an 
unsteered, 3x3 transmit array. 



AX 


AY 


1 H(f,r) 1 
( exact ) 


1 H(f,r) 1 
( simulated ) 


% error 


- 10 


10 


1.4324 e-2 


1.4318 e-2 


- 0.042 


20 


-20 


1.4324 e-2 


1.4320 e-2 


- 0.028 



Table 6.7 Acoustic field resulting from a 3 x 3 transmit array 
steered in the direction of the field point. 
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Naval Postgraduate School IBM 3033/4381 Network. It should be noted that the run 
times shown in Table 6.8 are for a frequency of 400 Hz. Note that the times detailed in 
Table 6.8 are for the simplest case possible. The phase integral term of the ocean medium 
transfer function has a closed-form expression for the isospeed medium case, which means 
that no numerical integration or approximation of the phase integral has to be made. 
Hence, computation time for this component of the integrand is at a minimum. The 
computation time becomes very significant when dealing with receive and / or transmit 
arrays rather than single elements. This is the main reason that most of our analysis has 
been carried out for small values of RLOS. 



RLOS 


Computation 


(m) 


Time ( min:sec ) 


10 


2:01 


100 


2:20 


1000 


3:23 


10000 


12:18 



Table 6.8 Simulation run times for a single point source 
to a single point receiver problem at 400 Hz. 
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VII. COMPUTER SIMULATION RESULTS FOR THE 
SINGLE GRADIENT MEDIUM CASE 



A. EVALUATION OF THE WKB PHASE INTEGRAL 

The results of the previous chapter validated the development carried out in Chapters 2 
through 4 for the case of a free-space, isospeed medium. The next step in the validation 
process is to apply the analysis to a more complex scenario, i.e., the free-space, single 
gradient case. For the purposes of this thesis we will characterize the sound-speed profile 
as being linear, with the speed of sound being given by 

c(y) = c yref + s( y ■ y ref)> (7-1) 



where g is the sound-speed gradient, and cyref 1S speed of sound at the depth y RE p 
The linear approximation of the sound-speed profile was selected as it allows for the easy 
description of a more complex medium by means of linear segments. The aim of this 
chapter is to validate the analysis for a single linear segment before the computer 
implementation is applied to a multiple linear segment sound-speed profile. 

Since the preceding development of the overall system complex frequency response 
given by Eq. (3.14) assumed a free-space problem, all we need to do is apply the condition 
of a single, linear, sound-speed profile to the ocean medium transfer function, given by 



H M (f,f r ,y 0 ;y) = 



K| k Y(y 0 )l p -J / y k Y(0 ^ e j27rf Y (y 0 Xy -y n ) 



2k Y (y 0 )V|k Y (y)| 



(4.34) 



126 



where 



k Y (y) = ±2*[(f/c(y)) 2 -f?] , fjfi(f/c(y)) 2 , 



(7.2a) 



and 




f r 2 >(f/c(y)) 2 . 



(7.2b) 



For the form of the sound-speed profile given by Eq. (7.1), we were unable to find a 
closed form solution to the phase integral in Eq. (4.34). This is in marked contrast to the 
results of the previous chapter, where the phase integral had a closed form expression 
which resulted in the ocean medium transfer function for the isospeed medium requiring no 
numerical integration in its evaluation. 

The method required to evaluate the phase integral of Eq. (4.34) is of critical 
importance. Recall from Chapter 3 that the expression for the overall system complex 
frequency response, given by Eq. (3.14), involved the evaluation of an integral whose 
integrand contained the ocean medium transfer function. In other words, there exists an 
integration within an integration. If a numerical integration routine, such as the IMSL10 
routine, DQDAG, discussed in the previous chapter, is used to evaluate this phase integral, 
then the resulting computational time involved is unmanageable. As an example, the 
computer simulation was run for a single frequency component at a RLOS of 100 m, and a 
single point source and a single point receiver. It was found for this case that the 
maximum allowable processing time of one hour on the Naval Postgraduate School IBM 
3033/4381 Network was attained before a solution could be computed. Consequently, we 
required an alternative method to numerical integration to evaluate the integral. 
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As an approximate solution to the phase integral, we expressed the propagation vector 
components in the Y (depth) direction, given by Eqs. (7.2a) and (7.2b), in terms of their 
binomial series. The reason for doing this was that the integral of each of the terms of the 
binomial series expansion, for both forms of the propagation vector component ky (y), has 
a closed form expression. As a result, the phase integral of the ocean medium transfer 
function can be replaced by a converging series. Taking the integral of the binomial series 
expansion of Eqs. (7.2a) and (7.2b), it can be shown that 



,,, , 27tf r i/frc(y)\ 2 i i/frc(y)\ 4 i 

,(y)dy = ± — |_ ta [c(y)] - ;r( — ) 2 '*(— ) 4 

i / f rC<y> \ 6 i s / f r c(y) \ 8 1 "I 
' 16 V f ) 6 ' 128 V f ) 8 ' J 



(7.3) 



for the propagating waves described by Eq. (7.2a), and 



f(y)dy + j2ic { - f r y - g [ 2 ( f r c(y) ) + 8 ( f r c(y) ) 3 
5 7 

+ 16 ( ?7cOO ) 5 + 128 ( fjx(y) ) 7 + ] } (7 4) 



for the evanescent waves described by Eq. (7.2b). The number of terms used in Eqs. 
(7.3) and (7.4) depends on the desired accuracy for the phase integral. Note that this 
method of determining the phase integral is inefficient for small values of ky (y). This is 
due to the fact that for small values of ky (y), f r » (f /c(y)) . As a result, the binomial 

expansion, and hence, the phase integral converge very slowly. 
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Recall from Chapter 3 that the plus (minus) sign in Eq. (7.3) is chosen whenever 
y > y Q (y < y G ) and corresponds to the field point being deeper (shallower) than the source 
point, i.e., downward (upward) traveling waves. The minus (plus) sign in Eq. (7.4) 
corresponds to the plus (minus) sign of Eq. (7.3) in order to generate evanescent waves. 
Recall that the full-wave solution can be thought of as an infinite summation of rays. If we 
think of rays from a single omnidirectional point source propagating in a medium 
characterized by a single, linear, sound-speed gradient, it can easily be visualized that some 
rays will never reach the receiver depth while others will reach the receiver depth twice. 
This situation is shown in Fig. 7.1 where it has been assumed that the receiver is below the 
source and the sound- speed gradient is positive. 

This is an important concept when we consider the way in which the overall system 
complex frequency response is calculated. Since the overall system complex frequency 
response is an integral with respect to fr and is a function of frequency, the WKB phase 
integral contained in the integrand is essentially the phase change of a particular ray, 
defined by the values of f and f r , traveling between the source and the receiver depth. For 
each of these rays, we need to decide how it behaves and determine a strategy for 
determining the value of the phase integral. 

This ray path problem was not encountered when dealing with the isospeed medium, 
since for this case the rays travel in straight lines as shown in Fig. 7.2. The situation 
shown in Fig. 7.2 assumes that the receiver is below the source and, as can be seen, all the 
downward traveling waves reach the receiver depth once in an isospeed medium. As a 
consequence, only the closed form expression of the phase integral for downward traveling 
waves needs to be used in the evaluation of the ocean medium transfer function for 
propagating waves. Similarly, only one form of the phase integral expression needs to be 
used for the evanescent waves. 
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4 




Fig. 7.1 Ray paths in a medium characterized by 
a single, positive, sound-speed gradient. 




( depth ) 

Fig 7.2 Ray paths in an isospeed medium. 
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Inspection of Fig. 7.1 shows the various groups of rays that will be encountered. 
Rays A and B are handled in the same manner. Since neither of these rays is refracted 
back to the receiver depth before reaching the horizontal range of the receiver, the phase 
integral is calculated for the downward traveling portion of the rays. For rays A and B, 
the integral is calculated from the source depth y Q meters to the first attainment of the 
receiver depth y meters, i.e., from y 0 to points 1 and 2, respectively. 

Ray C, unlike rays A and B, is refracted back to the receiver depth before the horizontal 
range of the receiver. As is illustrated in Fig. 7.1, the ray has both downward and upward 
portions that need to be taken into consideration. Consequently, the phase integral is 
calculated for a downward traveling ray from the source to the turning point at point 3, and 
for an upward traveling ray from the turning point back to the receiver depth at point 4. 
This essentially means that the phase integral has been subdivided into two separate 
segments, dependent on which form of the propagation vector given by Eq. (7.2a) is valid. 
Brekhovkikh [Ref. 6:pp. 66-67] also shows that such a ray encounters a the phase shift of 
-7t/2 radians at the turning depth, i.e., point 3. After passing through the turning point, the 
upward traveling ray keeps the -7i/2 phase shift. 

The final ray type to be considered is that characterized by ray D. As can be seen from 
Fig. 7.1, this propagating ray never reaches the receiver depth since it reaches a turning 
point above the receiver depth at point 5, and becomes an upward traveling ray. To 
evaluate the phase integral, we first consider the downward traveling portion of the ray 
down to the depth of the turning point. We then need to evaluate the integral from this 
point down to the receiver depth. There is clearly no propagating wave component going 
down to the depth of the receiver from the turning point, so any acoustic field reaching the 
depth of the receiver from this point must be evanescent. Accordingly, to evaluate the 
phase integral down to the receiver depth we consider the integrand to be the complex 
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propagation vector component given by Eq. (7.4). Recalling the discussion regarding the 
behavior of ray C at a turning point, this evanescent component can be thought of as what 
gives rise to the -tu/2 phase shift in the propagating portion of the ray reflected at the turning 
point. In other words, when the ray goes through a turning point, it gives rise to a phase 
shifted, reflected propagating component and a evanescent component 

Recall from the previous chapter the importance of evanescent waves. As was done 
for the isospeed medium case, we will use 



H(f,r) 



* / 0 



1.2(f/c 0 ) 



H M (f.fr.yo;y)Jo(2^rr m )e 



-j2jtf Y (y 0 Xy - y n ) 



frdfr. (6.5) 



as our approximate full-wave solution to the ocean propagation problem which means that 
we will be integrating with respect to fr up to a limit of 1.2 (f / c 0 ). Since this integration 
limit gave very good results for the isospeed medium case, it is assumed that it will also be 
applicable to the new single, linear gradient problem since the sound-speedprofile is slowly 
changing and does not exhibit any erratic behavior. Note that the region of integration 

(f/c 0 )<f r < 1.2(f/c 0 ) (7.4) 

describes evanescent waves resulting from the source condition and should not be confused 
with an evanescent wave that results from a ray being refracted through a turning point. 

B. RESULTS 

As discussed in the preceding section, the accuracy of the binomial series 
approximation depends on the number of terms used and for small values of ky(y) the 
number of terms required becomes very large. Just as for the situation in which we tried 
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to evaluate the phase integral using a numerical integration routine, if we specify a 
requirement that the phase integral be accurate to within 0.01% and evaluate the binomial 
expansion approximation for all the values of f r used to numerically integrate Eq. (6.5), the 
computation time becomes unmanageable. In this case it was again found that the 
maximum allowable processing time of one hour on the Naval Postgraduate School IBM 
3033/4381 Network was attained before a solution could be computed. 

As a consequence, it was decided to modify the region of integration in the evaluation 
of Eq. (6.5). Instead of integrating over the two regions 0 < f r < (f /c 0 ) and 
(f /c 0 ) < f r < 1.2 (f / c 0 ), we calculated the approximate full-wave solution of Eq. (6.5) 
using the following modified regions: 0 < f r < LIM (f / c 0 ) for initially propagating waves 
and (1/LIM) (f/c Q ) < f r < 1.2 (f/c 0 ) for initially evanescent waves. In these 
expressions, LEM is a constant less than unity that basically allows us to trade computation 
time against accuracy. After some investigation of the rate of convergence of the binomial 
series approximation of the phase integral, it was found that solutions could be obtained for 
a scenario in which the RLOS was 100 m, the transmitter was a single point source, and 
the receiver was a 3 x 3 array of point receivers, if the value of LIM did not exceed 0.995. 
It can easily be seen that the closer the value of LIM is to unity, the better our approximate 
full-wave solution should be. Accordingly, we conducted simulations for LIM equal to 
0.99 and 0.995. 

From our observations in the previous chapter, there are two forms of output that are 
useful in determining the performance of our full-wave approximation. First, there is the 
tabulation of the magnitude and phase of the complex acoustic field at the locations of the 
receive array elements. Recall that the usefulness of this form of output was that it allows 
for easy verification that the magnitude of the complex acoustic field is behaving correctly. 
For a given source / receiver geometry, it is an easy matter to determine if the the magnitude 
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of the complex acoustic field is increasing or decreasing according to theory as we go from 
one element to the next. Second, there is the output from the adaptive beamforming 
algorithm, used in the previous chapter, processing the time-domain electrical signal at each 
element of the receive array. This algorithm primarily uses the phase of the received signal 
in its processing. Therefore, it validates how accurate the phase of the predicted complex 
acoustic field is. Consequently, analysis of both these outputs provides a check of the 
accuracy of both the amplitude and phase of the complex acoustic field at a specified field 
point. 

Computer simulation results were generated for two values of the constant LIM, as 
previously discussed, in conjunction with a source / receiver geometry of the form shown 
in Fig. 7.1. The source was placed at a depth of 1,000 m, the receiver was below at a 
depth of 1,050 m, and the RLOS was 100 m. The sound-speed profile of the ocean 
medium was characterized by the linear sound- speed profile given by Eq. (7.1), where 
yREF = 1.000 m, cyref = 1.475 m/sec and the gradient had the value 0.017 1/sec. The 
aim of investigating these particular values of LIM was to determine whether the small 
region of integration that was discarded was significant. 

The values of the magnitude and phase (in degrees) of the acoustic field for four test 
cases are shown in Tables 7.1 through 7.4. Note that the first value is the magnitude and 
the second value is the phase. Given the placement of the receive array relative to the 
source, i.e., XR = 0 and YR = 1050 m, which is below the source, we would expect the 
magnitude of the acoustic field to decrease with increasing n. This corresponds to reading 
down the table or equivalently increasing the depth of the receive element or field point. 
We would also expect the value of the acoustic field to decrease with increasing absolute 
value of m. This corresponds to reading across the table out from the center or 
equivalently increasing the displacement in the X direction of the element or field point. 
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We expect the magnitude of the acoustic field to decrease in both of these situations since 
their effect is to increase the radial range between the source and the receive element , or 
field point Note that these are the same results as expected for the isospeed medium case. 
This is because the RLOS is a small value, and hence, the propagating rays have not had 
much of an opportunity to refract Therefore, a crude approximation is that the medium is 
homogeneous over this short range. 

Case A, shown in Table 7.1, evaluates the overall system complex frequency response 
over the modified region of the initially propagating waves. In other words, we only 
integrate over the region 0 < f r < UM (f / c 0 ) and ignore the region of initially evanescent 
waves. Case B, on the other hand, shown in Table 7.2, is the same as case A except that 
it also includes the modified region of intitially evanescent waves. Both cases are 
computed for LIM = 0.99. It can be easily observed from these two tables that the 
magnitude of the acoustic field does not behave exactly as expected. In both cases, the 
magnitude of the acoustic field decrease as we read down the table as expected. However, 
in both cases, the second and third rows have the magnitude increasing as we move away 
from m = 0. 

Note that the magnitude and phase of the complex acoustic field is very similar for both 
these cases despite the inclusion of evanescent waves in Case B. At first, this may appear 
to contradict the results of the previous chapter. For the isospeed medium the inclusion of 
evanescent waves had a dramatic effect on the value of the complex acoustic field. So 
what causes this apparently contradictory behavior? The difference between the observed 
results for the isospeed and single sound-speed gradient cases is due to the region of 
evanescent waves considered. The evanescent waves in the region very close to 
f r = (f /c 0 ) are by far the most dominant ones. Accordingly, the modified region of 
evanescent waves ignores the most dominant ones. So why the upper limit of 
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n/m 


1 


0 


-1 




8.17543 e-4 
- 165.452 


8.18198 e-4 
- 164.725 


8.17543 e-4 
- 165.452 


0 


8.12795 e-4 
139.886 


8.12545 e-4 
140.616 


8.12795 e-4 
139.886 


1 


7.98127 e-4 
83.940 


7.97160 e-4 
84.627 


7.98127 e-4 
83.940 



Table 7.1 Case A : LIM = 0.99, initially evanescent waves not included. 



n/m 


1 


0 


-1 


-1 


8.17524 e-4 
- 165.451 


8.18179 e-4 
- 164.725 


8.17524 e-4 
- 165.451 


0 


8.12783 e-4 
139.885 


8.12534 e-4 
140.615 


8.12783 e-4 
139.885 


1 


7.98128 e-4 
83.939 


7.97161 e-4 
84.626 


7.98128 e-4 
83.939 



Table 7.2 Case B : LIM = 0.99, initially evanescent waves included. 
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f r = 1.2 (f / c 0 )? This upper limit becomes significant when considering very short 
ranges, i.e., ranges less than 100 m. At these short ranges the region of significant 
evanescent w-aves increases rapidly with decreasing range. The specific upper limi t of 
f r = 1.2 (f /c 0 ) w as chosen to give accurate values for the acoustic field calculated from 
the approximate full-w’ave solution down to a range of 10 m. 

If our reasoning has been correct, then we would expect our values for the computed 
complex acoustic field to improve as we increased the value of LIM closer to unity. Case 
G, which ignores the modified region of evanescent w-aves, and Case D, which includes 
the modified region of evanescent waves, were both computed for LIM = 0.995 and are 
shown in Tables 7.3 and 7.4, respectively. As can be seen from the values in these tables, 
the same phenomena is demonstrated as was evident when comparing Tables 7.1 and 7.2. 
In fact, it can be easily seen that increasing the region of evanescent w'aves has had minimal 
effect on the value of the complex acoustic field. Does this mean that our explanation of 
the importance of the evanescent waves in the region very close to f r = (f/c 0 ) is 
incorrect? Let us reconsider the isospeed medium of the previous chapter and see what 
happens w hen w e apply these modified regions of integration. 

Consider the same source / receiver geometry as we have been using for this 
investigation, but apply the approximate full- wave solution to an isospeed problem, where 
the constant speed of sound c© is 1,475 m/sec as was the case in the previous chapter. 
Applying the modified regions of integration, with LIM = 0.999, to the ocean medium 
transfer function given by Eq. (6.5) results in the values of the complex acoustic field 
shown in Table 7.5. Note the way that the magnitudes of the acoustic field do not behave 
as predicted In the third row the magnitude of the acoustic field increases as w e move 
away from m = 0. Also, note that the magnitudes of the acoustic field increase with 
increasing depth. Both of these results are the opposite from what we expect. Comparing 
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n/m 


1 


0 


-1 


-* 


8.29299 e-4 
- 162.942 


8.28427 e-4 
- 162.348 


8.29299 e-4 
- 162.942 


0 


7.94735 e-4 
142.244 


7.94784 e-4 
142.804 


7.94735 e-4 
142.244 


1 


7.61585 e-4 
84.539 


7.62560 e-4 
85.126 


7.61585 e-4 
84.539 



Table 7.3 Case C : LIM = 0.995, initially evanescent waves not included. 



n/m 


1 


0 


•- 


-1 


8.29110 e-4 
- 162.938 


8.28241 e-4 
- 162.344 


8.29110 e-4 
- 162.938 


0 


7.94604 e-4 
142.237 


7.94653 e-4 
142.797 


7.94604 e-4 
142.237 


1 


7.61599 e-4 
84.529 


7.62572 e-4 
85.116 


7.61599 e-4 
84.529 



Table 7.4 Case D : LIM = 0.995, initially evanescent waves included. 
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Tables 7.5 and 6.5, we observe that there is a significant difference between the two sets of 
results, yet the only difference between the two was that a modified region of integration 
was used to calculate the values in Table 7.5. Recall that the value of LIM for the values 
calculated in Table 7.5 was 0.999, that is, very close to unity. These results reaffirm our 
earlier conclusions that the evanescent waves in the region very close to f r = (f/c 0 ) 
contribute very significantly to the overall value of the integral. This implies that around 
this value of f r> there may be a region of relatively constant phase in the integrand of the 
overall system complex frequency response. 

The outputs from the LMS adaptive beamforming algorithm for Cases A through D are 
shown in Table 7.6. Recall the geometry illustrated in Fig. 6.3. The angle definitions are 
still valid with the exception that they relate to the local angles of arrival of the wave at the 
receive transducer. In other words, for the isospeed medium case the local angle of arrival 
was also the direction of the source since all rays traveled in straight lines. From the 
values for all four cases we see that the range estimates RLOS (X) are in error by an order 
of 10 - 15 % while the estimates RLOS (Y) are in error by up to 65 %. This phenomena 
was also observed for the isospeed case when the effects of evanescent were ignored. For 
that case, it was argued that since the overall system complex frequency response is 
axisymmetric, then the relative behavior of the acoustic field in the X direction should be 
minimally affected by the behavior of the acoustic field in the Y direction. Consequently, 
if there is a problem with the description of the full-wave solution in the Y direction, as in 
this case, the value of RLOS (X) should still be fairly accurate, which is the result that we 
observed. 

The angular outputs shown in Table 7.6 agree very closely with theory. Even though 
the angles are the predicted local angle of arrival, for our case where the RLOS is small 
(i.e., 100 m), the effects of refraction will have been small and the local angle of 
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n/m 


1 


0 


-1 




7.33866 e-4 
- 169.589 


7.34792 e-4 
- 168.589 


7.33866 e-4 
- 169.321 


0 


7.75856 e-4 
131.862 


7.77764 e-4 
132.604 


7.75856 e-4 
131.862 


1 


8.30442 e-4 
76.103 


8.29309 e-4 
76.781 


8.30442 e-4 
76.103 



Table 7.5 Case E : LIM = 0.999, initially evanescent included. 



CASE 


LIM 


RLOS (X) 
(m) 
est 


RLOS (Y) 
(m) 
est 


© 

(deg) 

est 


(deg) 

est 


A 


0.99 


90.690 


75.577 


29.471 


270.00 


B 


0.99 


90.688 


75.628 


29.471 


270.00 


C 


0.995 


115.748 


34.343 


30.015 


270.00 


D 


0.995 


115.720 


34.438 


30.019 


270.00 



Table 7.6 Output from adaptive LMS beamforming algorithm. 
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arrival should be very close to the angles describing the position of the source relative to the 
receiver. Any effect due to refraction for the given geometry and sound-speed profile will 
be to decrease the angle 0 slightly and leave 'F unaffected. However, our results are too 
inconsistent to tell if there is a noticeable decrease in the value of 0, but they do 
demonstrate that remains unchanged. 

The results obtained from this investigation of a medium characterized by a linear 
sound-speed profile with a single gradient have essentially been inconclusive. It was 
found that the evaluation of the phase integral in the ocean medium transfer function based 
on the WKB approximation is of crucial importance to the performance of the approximate 
full-wave solution. If the linear systems approach, and the coupling equations developed 
in this thesis are to be successfully applied to the WKB approximation of the ocean 
medium, we require a fast, accurate method of determining the value of the phase integral 
to be able to accurately compute the value of the complex acoustic field at some given 
spatial location. 

The method of representing the phase integral in terms of a power series expansion, 
while better that using numerical integration techniques, was inappropriate since there is a 
major contribution to the value of the acoustic field in the region where this expansion 
converges very slowly. The use of the series expansion allowed approximate values of the 
acoustic field to be calculated, whereas the use of the IMSL10 numerical integration routine 
DQDAG took too much computing time and no useful results could be obtained. 
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vm. CONCLUSIONS AND RECOMMENDATIONS 



It appears from our preliminary results that our approach to modelling the ocean 
medium pulse propagation problem, which incorporates both linear systems theory and the 
physics of acoustic propagation, is a valid one. 

A signal generation scheme based on a truncated Fourier series and complex envelope 
representation of narrowband signals was implemented and shown to work correctly. The 
scheme is capable of generating CW and LFM pulses with rectangular, Hamming or 
Hanning amplitude weighting suitable for conversion into a transmitted acoustic signal. 

We have shown that the beam pattern of the transmit array does not violate axial 
symmetry. This is a result of the theory of linear superposition, which allows us in the 
case of our linear system model, to decompose the beam pattern into a sum of 
omnidirectional point sources. 

The model has been extensively tested and validated for the case of a free-space 
isospeed medium. It was found that the formal solution to the propagation problem could 
be modified to include only a small portion of the region of evanescent waves. This 
approximate full-wave solution has demonstrated that it gives extremely accurate results for 
radial ranges between source and receiver down to 10 m. It was also shown that the 
formal solution was inaccurate in the situation where the source and field point were at or 
very near the same depth. It was concluded that this was due to the nature of the WKB 
approximation rather than a problem with the system model. It is recommended that a 
modified WKB approximation be used in the vicinity of a turning point, which involves 
Airy functions. This modified solution should be investigated in an attempt to determine if 
it improves the model performance for this source / receiver geometry. 
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Inconclusive results were obtained about the model’s behavior for an ocean medium 
characterized by a linear sound-speed profile with a single gradient. This was due to the 
difficulties encountered in evaluating the phase integral term of the ocean medium transfer 
function based on the WKB approximation. For this case we were unable to determine a 
closed form expression for the phase integral. Numerical integration techniques were too 
slow and our derived approximation to the integral was unsuited to the numerical behavior 
of the integrand and its effect on the behavior of the approximate full-wave solution to the 
propagation problem. Further investigation needs to be carried out on how to quickly and 
accurately evaluate this integral. Some recommended avenues for further investigation 
are : 



• to determine whether a different description of the sound-speed profile will yield a 
closed form expression for the phase integral, 

• to find a better approximation of the phase integral than the binomial series 
expansion. 

There is also a question regarding the strategy required to implement the phase integral. 
From a ray acoustics viewpoint, this integral represents the phase change experienced by a 
ray as it travels from the source to the receiver depth. For the inhomogeneous medium, 
these rays will be refracted and there exists the situation were the ray doesn't reach the 
receiver depth at all, or reaches it more than once. How should the integral be evaluated in 
these circumstances? Because of the inconclusive nature of our computer simulation 
results, we were unable to validate or disprove our scheme to compute the phase integral. 

It was demonstrated that the outputs from the computer simulation, i.e., the magnitude 
and phase of the complex acoustic field as a function of frequency and spatial coordinates, 
and the time-domain output electrical signal from each element in a receive array were 
implemented correctly. While no matched-field processing was done using this first 
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output, it was demonstrated that the plotting of the magnitude of the complex acoustic field 
was a powerful aid in helping to visualize what the acoustic field looked like. The time- 
domain signals were processed by an external frequency domain adaptive beamforming 
program and the data obtained was very useful in validating the performance of the model. 
This also demonstrated that generated data could be successfully passed to an external file, 
to be processed by an ancillary program of our choosing. Thus, the generated data is 
independent of whatever receive array signal processing we wish to carry out 

It has been shown that an approximate solution to the linear wave equation can be 
implemented as an ocean medium transfer function. However, except when closed form 
solutions or short recursive computation schemes exist for all terms in the ocean medium 
transfer function, we are constrained by the required computation time. It is recommended 
that the feasibility of the following be investigated to decrease computation time : 

• optimization of the computer code at the expense of program modularity, 

• parallel processing techniques, 

• availability of a faster computing system than the IBM 3033/4381 Network. 

The very nature of the full-wave solution to the pulse propagation problem makes our 
model computation intensive. If this method of solution is to be a useful tool in 
investigating realistic scenarios, we need to consider either increasing the speed with which 
the current calculations are handled, or use another technique, such as the Fast-Field- 
Program method of evaluating the wave vector component integrals. 
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